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FOREWORD 


A  pilot-induced  oscillation  (PIO)  results  from  the  interaction  of  the  pilot  and  the  dynamics 
of  the  vehicle  being  controlled.  It  may  be  caused  or  affected  by  several  elements  of  the 
aircraft  design  or  the  mission  task.  PIO  affects  the  pilot’s  ability  to  perform  a  given  task, 
ranging  from  an  annoying  aircraft  motion  to  inability  to  complete  the  task  to,  in  the  most 
extreme  cases,  jeopardizing  the  safety  of  the  aircraft  and  crew.  Because  it  occurs 
sporadically,  PIO  can  be  one  of  the  most  insidious  flying  qualities  problems. 

Criteria  currently  exist  for  some  of  the  elements  of  aircraft  design  which  influence  PIO 
tendencies;  however,  the  criteria  in  the  flying  qualities  standard  (MIL-STD-1797)  address 
the  effect  of  these  elements  individually  rather  than  cumulatively.  These  elements  include 
sluggish  response  modes,  low  damping,  excessive  phase  lag  or  time  delay,  overly  sensitive 
stick  gradients,  unstable  response  modes,  coupled  response  modes,  etc.  Furthermore, 
there  are  some  elements  that  are  known  to  influence  PIO  tendencies  which  are  not 
addressed  in  MIL-STD-1797.  Some  examples  of  these  elements  are  aerodynamic 
nonlinearities,  control  system  nonlinearities,  actuator  rate  limits,  nonlinear  stick  gradients, 
etc. 

This  report  is  the  first  step  of  an  Air  Force  initiative  to  develop  a  comprehensive  and 
unifying  theory  to  explain  the  nature  of  the  interaction  of  pilot  and  vehicle  which  results  in 
PIO  and  to  develop  the  capability  to  predict  an  aircraft’s  PIO  tendencies  due  to  the 
combined  effect  of  all  influencing  elements.  This  theory  must  be  consistent  with  several 
PIOs  which  have  been  documented  in  research,  development,  and  operational  aircraft. 

This  effort  should  also  produce  design  criteria  which  can  be  used  in  the  development 
process  to  assess  the  risk  of  PIO  prior  to  manned  simulation. 

The  four  volumes  of  this  report  represent  a  broad-band  approach  to  developing  this 
comprehensive  theory.  There  were  five  contractors  involved:  Calspan  Corporation;  Hoh 
Aeronautics,  Inc.;  McDonnell  Douglas  Aerospace;  Systems  Technology,  Inc.;  and  Virginia 
Polytechnic  Institute  and  State  University.  The  objectives  of  the  contractors  were  to 
develop  their  theories  and  partially  validate  them  against  existing  data. 

Volume  I  is  the  work  of  Systems  Technology,  Inc.  This  research  included:  compilation  of 
available  PIO  time  histories  and  references  as  an  initial  step  toward  a  comprehensive  PIO 
database;  refinement  of  the  proposed  PIO  categories  defined  in  NASA  CR-4683; 
development  of  PIO  theories  based  on  these  categories;  development  of  methods  to 
handle  the  higher  frequency  and  nonlinear  aspects  of  PIO  analysis  with  an  emphasis  on 
rate  limiting;  and  a  review  of  existing  and  proposed  linear  PIO  criteria. 

Volume  n  is  the  combined  work  of  McDonnell  Douglas  Aerospace,  Advanced  Transport 
Aircraft  at  Long  Beach,  CA,  McDonnell  Douglas  Aerospace,  Advanced  Programs,  at  St. 
Louis,  MO,  and  Hoh  Aeronautics,  Inc.  at  Lomita,  CA.  There  are  primarily  two 
components  associated  with  this  study.  One  element  is  an  examination  of  PIO  events  that 
have  occurred  in  the  course  of  initial  flight  development  on  several  aircraft  produced  by 


the  McDonnell  Douglas  Corporation.  The  other  element  is  an  exploration  of  several 
theories  associated  with  PIOs,  particularly  those  with  rate  limiting  involved. 

Volume  in  was  produced  by  the  Virginia  Tech  Department  of  Aerospace  and  Ocean 
Engineering.  This  report  describes  an  analysis  mediod  capable  of  predicting  PIO 
tendencies  due  to  several  simultaneous  factors.  A  unified  approach  involving  pilot 
modeling,  stability  robustness  analysis,  and  multivariable  describing  function  analysis  is 
applied  to  the  problem  of  identifying  aircraft  with  PIO  tendencies.  The  approach  is  shown 
to  have  ties  to  existing  PIO  criteria  and  is  successfully  applied  to  the  prediction  of  PIO 
tendencies  of  the  M2-F2  lifting  body. 

Volume  rV  is  the  result  of  Calspan’s  efforts.  This  report  presents  the  theory,  fundamental 
principles,  and  analytical  procedures  of  a  quantitative  criterion  for  the  prediction  of  PIO 
tendencies  based  on  a  “time-domain  Neal-Smith”  approach.  The  criterion  is  validated 
against  three  very  reliable  flying  qualities  data  bases.  At  present  the  criterion  is 
intentionally  limited  to  the  evaluation  of  pitch  control  only.  No  fundamental  limitations 
were  discovered  which  preclude  the  evolution  of  this  mediodology  and  analytical 
procedures  to  PIO  analysis  of  roll  control  or  “outer-loop”  longitudinal  control,  such  as 
control  of  aircraft  flight  path. 

After  further  evaluation  and  deliberation  of  these  reports,  the  Air  Force  will  pursue  the 
most  promising  ideas,  or  combination  of  ideas,  and  validate  them  with  further  simulation 
and  flight  testing. 
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PREFACE 


A  pilot-induced  or  pilot-in-the-loop  oscillation  (PIO)  is  a  very  complicated 
phenomenon  stemming  from  a  dynamic  interaction  between  the  pilot  and  the  aircraft.  This 
interaction  means  that  the  aircraft  can  be  otherwise  well  behaved  and  the  pilot  can  be  very 
well  trained;  however,  when  tight  interaction  between  the  pilot  and  the  aircraft  is  required, 
an  undesired  oscillation  can  result.  Furthermore,  the  adaptive  nature  of  the  human  pilot 
makes  such  oscillations  very  difficult  to  predict.  This  report  describes  an  analysis  method 
capable  of  predicting  PIO  tendencies  due  to  several  simultaneous  dynamic  factors.  A 
unified  approach  involving  pilot  modeling,  stability  robustness  analysis,  and  multivariable 
describing  function  analysis  is  applied  to  the  problem  of  identifying  aircraft  with  PIO 
tendencies.  The  approach  is  shown  to  have  ties  to  existing  PIO  criteria  and  is  successfully 
applied  to  the  prediction  of  PIO  tendencies  of  the  M2-F2  lifting  body.. 
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1  Introduction 


A  pilot-induced  or  pilot-in-the-loop  oscillation  (PIO)  is  a  very  complicated  phenomenon 
stemming  from  a  dynamic  interaction  between  the  pilot  and  the  aircraft.  This  interaction  means 
that  the  aircraft  can  be  otherwise  well  behaved  and  the  pilot  can  be  very  well  trained;  however, 
when  tight  interaction  between  the  pilot  and  aircraft  is  required,  an  undesired  oscillation  can 
result.  Furthermore,  the  adaptive  nature  of  the  human  pilot  makes  such  oscillations  very  difficult 
to  predict. 

One  of  the  most  recent  aircraft  loses  due  to  a  PIO  was  the  Lockheed/Boeing  YF-22A 
advanced  tactical  fighter  prototype.^]  Figure  1-1  shows  a  re-creation  of  the  PIO  that  occurred  in 
the  YF-22A.  This  figure  was  created  electronically  scanning  the  last  seven  seconds  of  the  time 
history  of  Reference  [1].  While  the  figure  is  not  drawn  completely  to  scale  (the  altitude 
excursions  have  been  accentuated),  the  figure  does  illustrate  the  rather  dramatic  changes  in  flight 
attitude  that  can  occur  during  a  fiilly  developed  PIO. 

The  YF-22A  prototype  crash  is  an  example  of  a  very  complicated  dynamic  situation  that  led 
to  a  PIO.  Several  significant  changes  occurred  simultaneously  in  the  aircraft  dynamics. 

According  to  Reference  [I],  the  landing  gear  of  the  aircraft  were  retracted  at  nearly  the  same  time 


Figure  1-1.  YF-22A  Pilot-Induced  Oscillation. 
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as  the  afterburner  was  ignited  to  initiate  a  planned  go-around.  These  configuration  changes 
caused  significant  mode  switching  in  the  aircraft  control  laws.  The  control  law  changes  are 
usually  governed  by  "transient-free"  switches  that  slow  the  transition  of  control  system  gains  from 
landing  to  up-and-away  configuration.  The  pilot,  on  the  other  hand,  must  adapt  his  control 
strategy  continuously  and  at  the  same  rate  as  the  aircraft  dynamics  are  changing.  Therefore,  the 
PIO  of  the  YF-22A  is  an  example  of  several  different  dynamic  elements  contributing  to  the 

oscillation. 

The  increased  complexity  of  modern  aircraft  appears  to  increase  the  likelihood  of  PIO  due  to 
several  interacting  dynamic  factors.  Since  the  YF-22A  prototype  was  highly  instrumented  and  the 
PIO  event  was  captured  on  film,  it  was  possible  to  deduce  the  various  dynamic  phenomenon  that 
ultimately  lead  to  the  crash. 

Could  any  of  the  existing  PIO  analysis  methods  have  predicted  the  YF-22A  PIO  before  it 
occurred?  Most  existing  PIO  analysis  methods  and  criteria  were  developed  to  address  individual 
dynamic  effects.  For  example,  it  is  usually  assumed  that  the  pilot  is  only  controlling  a  single  axis 
or  moving  only  one  manipulator.  Other  criteria  cannot  handle  nonlinearities  and  few,  if  any,  can 
address  the  combined  effect  of  several  dynamic  factors  as  in  the  YF-22A  crash.  Consequently,  it 
is  unlikely  that  an  existing  analysis  method  could  have  predicted  the  YF-22A  PIO  incident. 

1.1  Scope 

This  report  deals  primarily  with  computational  analysis  techniques  designed  to  predict  PIO 
susceptibility  of  modern  aircraft.  These  techniques  are  intended  for  use  during  aircraft 
development  or  modification  prior  to  flight.  Methods  to  study  PIO  behavior  during  flight  test  and 
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manned  simulations  will  not  be  discussed.  Aircraft  and  pilot  modeling  will  only  be  addressed  to 
the  extent  required  to  utilize  the  PIO  analysis  methods. 

1.2  Objective 

A  comprehensive  and  unified  PIO  analysis  method  is  sought.  This  unified  method  must  be 
able  to  handle  the  combined  effect  of  all  influencing  dynamic  elements.  The  theory  must  also 
explain  the  nature  of  interaction  of  the  pilot  and  vehicle  that  causes  the  PIO.  And  finally,  it  must 
be  consistent  with  several  PIOs  which  have  been  documented  over  the  last  few  years. 

The  basic  premise  of  this  research  is  that  a  PlO-resistant  aircraft  is  resistant  to  closed-loop 
pilot  dynamic  variations  in  the  face  of  nonlinear  aircraft  behavior.  It  could  be  argued  that 
nonlinear  aircraft  behavior  causes  the  pilot  to  change  his/her  equalization  instead  of  vice  versa 
From  an  analysis  point-of-view,  the  simultaneous  action  of  the  aircraft  and  pilot  dynamics  are 
addressed  from  either  viewpoint. 

The  possible  variations  in  pilot  dynamics  that  might  lead  to  or  "trigger"  a  PIO  include 
primarily  an  increase  in  pilot  gain  due  to  an  excited  or  overly  aggressive  pilot.  In  addition, 
aggressiveness  may  also  alter  the  pilot's  inherent  time  delay  or  phase  lead  generating  capability. 
Human  operators  are  also  known  to  exhibit  nonlinear  behavior,  such  as  "bang-bang"  control,  l^l 
Any  unified  PIO  analysis  method  should  be  able  to  address  any  or  all  of  these  variations. 

The  principal  nonlinear  aircraft  dynamic  effect  that  is  known  to  contribute  to  PIO  is  control 
surface  actuator  rate  and  deflection  limiting.  Other  effects  such  as  control  system  mode  switching 
and  nonlinear  stick  gradients  are  known  to  contribute  to  PIO  susceptibility.  Some  PIO  events 
have  lead  to  such  large  amplitude  oscillations  that  nonlinear  aerodynamics  have  played  a 
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significant  role.  A  unified  PIO  analysis  method  should  be  able  to  handle  each  of  these  dynamic 
elements  as  well  as  combinations  of  them. 

1.3  Plan 

Since  the  basic  premise  of  this  research  involves  the  study  of  pilot/aircraft  stability  robustness 
and  nonlinear  system  oscillations,  the  unified  analysis  method  proposed  is  a  combination  of 
multivariable  stability  robustness  analysis  and  multivariable  describing  function  analysis.  These 
methods  were  developed  specifically  to  handle  multiple  dynamic  effects  simultaneously.  In 
addition  to  this  final  report,  the  following  work  plan  was  devised. 

The  first  task  consists  of  the  development  of  the  necessary  stability  analysis  software. 

Existing  computational  analysis  algorithms  have  been  collected  and  refined.  The  MATLAB^“ 
computer-aided  engineering  software  package  is  used  as  the  computational  environment  in  which 
all  analysis  is  conducted.  Since,  other  competing  packages  have  similar  programming  (or  macro) 
features,  it  is  expected  that  translation  to  other  packages  can  be  easily  completed. 

The  four  basic  computational  tasks  that  must  be  carried  out  are:  1)  compute  the  pilot  model, 
2)  construct  the  block  diagram,  3)  perform  a  structured  singular  value  analysis,  and  4)  search  for 
limit  cycles.  Algorithms  already  exist  (in  MATLAB™  form)  to  compute  optimal  control  models 
of  the  human  pilot  [3],  and  no  modification  is  necessary  for  use  in  the  current  PIO  analysis. 
SIMULINK^^,  an  optional  program  linked  to  MATLAB™,  is  used  to  perform  the  necessary  block 
diagram  manipulations  and  simulations.  SIMULINK™  is  an  excellent  environment  for  study  of 
systems  with  simple  isolated  nonlinearities  of  the  type  generally  used  in  PIO  investigations.  The 
algorithms  to  perform  the  real  structured  singular  value  analysis  and  the  limit  cycle  search  were 
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also  coded  within  the  computational  framework  of  MATLAB™.  The  algorithms  are  based  on 
those  currently  available  in  the  open  literature. 

The  second  task  was  to  perform  verification  case  studies.  Once  all  of  the  analysis  codes  have 
been  implemented  and  tested,  case  studies  are  conducted  using  existing  aircraft  data.  The  HAVE 
PIO  database  is  used  to  study  existing  criteria  using  multivariable  stability  robustness  methods. 
This  activity  is  intended  to  show  that  the  new  unified  PIO  analysis  methods  can  be  readily  tied  to 
existing  methods  using  an  existing  database. 

The  multivariable  limit  cycle  software  is  separately  validated  using  a  previous  study  of  back¬ 
up  control  system  designs  for  the  F-4  aircraft.  The  original  study  was  intended  to  show  the  PIO 
susceptibility  of  the  back-up  control  systems.  This  application  is  particularly  appealing  because 
both  control  surface  rate  and  deflection  limiting  were  shown  to  occur  simultaneously. 

The  third  task  was  to  develop  new  design  criteria  from  the  new  PIO  theory.  The  HAVE  PIO 
database  is  used  again  to  develop  a  new  single-loop  PIO  criteria  from  multivariable  stability 
robustness  theories.  The  M2-F2  aircraft  is  used  to  study  the  effect  of  combinations  of  nonlinear 
dynamic  elements.  The  M2-F2  application  illustrates  how  the  unified  PIO  analysis  method  can  be 
used  to  isolate  combinations  of  dynamic  elements  which  lead  to  oscillatory  behavior. 
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2  Background 


Aircraft  characteristics  that  are  known  to  contribute  to  pilot-in-the-loop  oscillation  tendencies 
include;  sluggish  response  modes,  low  damped  modes,  excessive  phase  lag  or  time  delay, 
sensitive  stick  gradients,  unstable  response  modes  (that  cannot  be  stabilized  by  the  pilot),  and 
unusual  coupling  responses.  These  characteristics  have  been  studied  by  researchers  since  the 
early  1960's.  The  existing  research  literature  and  the  current  military  flying  qualities  speciflcation 
present  analysis  methods  that  have  been  successful  in  predicting  pilot-in-the-loop  oscillation 
behavior  for  each  of  these  characteristics.  However,  with  the  exception  of  Reference  [4],  very 
little  research  has  been  devoted  to  the  study  of  several  of  these  dynamic  effects  occurring 
simultaneously. 

Early  research  concerning  PIO  behavior  is  typified  by  the  paper  by  Hirsch  and  McCormick,  f^l 
Hirsch  and  McCormick  used  fixed-  and  motion-based  simulations  to  study  generalized 
longitudinal  aircraft  configurations.  However,  their  inclusion  of  simple  pilot  compensation 
strategies  was  recognized  early.  Their  analysis  consisted  of  experimental  observations  and  simple 
root-locus  charts.  In  addition,  Hirsch  and  McCormick  revealed  the  criticality  of  including  both 
visual  and  motion  cues  in  their  analysis. 

The  PIO  theory  developed  by  R.  Smith  is  perhaps  the  most  widely  recognized  analysis  and 
prediction  method.  1^1  Smith's  method  is  described  in  the  military  flying  qualities  specification  and 
has  been  evaluated  by  other  researchers,  In  the  authors'  view,  the  success  of  Smith's  method 
is  due  to  the  fact  the  interaction  of  the  pilot  and  aircraft  is  properly  modeled  using  elements  of 
McRuer's  crossover  pilot  model.  I '•’1  Thus,  the  criterion  is  based  not  just  on  the  dynamics  of  the 
aircraft,  but  also  on  the  human  pilot  capabilities. 
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One  of  the  most  prolific  researchers  in  the  area  of  PIO  investigation  has  been  R.A.  Hess.[*‘-i-^l 
Dr,  Hess  has  used  much  more  complex  and  intricate  models  of  the  human  operator  in  an  effort  to 
more  completely  represent  the  interaction  between  the  pilot  and  the  aircraft.  Hess's  1991  paper. 
Reference  [11],  develops  a  PIO  prediction  method  that  limits  the  bandwidth  and  slope  of  the 
pilot/vehicle  open-loop  transfer  fiinction.  The  success  of  Hess's  method  therefore  relies  heavily 
on  an  accurate  model  of  the  human  pilot. 

In  his  PIO  studies,  Hess  has  been  using  the  optimal  control  model  (OCM)  of  the  human 
operator.  This  model  was  developed  by  KJeinman,  Baron,  and  Levison  in  1970.n41  It  has  also 
been  used  by  other  researchers  to  develop  PIO  criteria  for  different  aircraft  types  and  flight 
tasks.l‘5-'7]  Like  Hess's  criteria,  these  researchers  have  developed  criteria  using  frequency- 
domain  representations  of  either  the  open-  or  closed-loop  pilot  vehicle  transfer  functions.  Closed- 
loop  resonant  peaking  behavior  (References  [15]  and  [16])  and  open-loop  crossover  properties 
(References  [13]  and  [17])  have  been  effective  PIO  predictors. 

Some  research  has  been  devoted  to  design  methods  to  reduce  PIO  tendencies.  However, 
since  most  of  the  prediction  methods  have  close  connections  to  feedback  control  theory,  aircraft 
control  system  remedies  are  usually  evident.  In  cases  where  flight  control  system  modifications 
are  unwise  or  costly.  References  [18],  [19],  and  [20]  describe  ways  in  that  the  control  stick 
dynamics  can  be  altered  to  help  alleviate  or  suppress  PIO  behavior. 


7 


3  PIO  Theory 


Before  the  specifics  of  the  proposed  analysis  method  are  set  forth,  some  basic  concepts 
concerning  PIO  are  addressed.  The  working  definition  of  a  PIO  used  in  this  report  is  given.  The 
acceptable  and  unacceptable  characteristics  of  PIO  are  then  stated.  Finally,  some  of  the  known 
and  suspected  causes  of  PIO  are  listed. 

3.1  Definitions 

Much  discussion  has  been  and  probably  will  continue  to  be  fervently  pursued  on  the  precise 

definition  of  a  PIO.  For  the  purposes  of  this  study,  a  simple  definition  is  employed, 

A  PIO  results  when  the  closed-loop  pilot/vehicle  system  oscillates  due  to  a  loss  in  asymptotic 
stability. 

Since  the  current  analysis  methods  deal  primarily  with  mathematical  models  of  the  pilot  and 
aircraft,  this  definition  is  easily  tested.  Asymptotic  stability  of  a  linear  system  is  lost  when  any 
pole  (or  complex  pair  of  poles)  cross  from  the  left  half  plane  onto  the  jco  axis.  For  a  nonlinear 
system,  in  addition  to  examining  the  poles  of  the  corresponding  linearized  system,  the  existence  of 
a  limit  cycle  (sustained  oscillation)  will  be  taken  as  a  separate  measure  that  precludes  the 
condition  of  asymptotic  stability  of  the  closed-loop  system. 

3.2  Acceptable/Unacceptable  Characteristics 

The  term  PIO  is  certainly  recognized  by  all  those  who  work  in  the  design  and  flight  testing  of 
high  performance  military  as  well  as  civil  transportation  aircraft.  If  we  proceed  from  the 
standpoint  that  a  PIO  prone  aircraft  is  undesirable,  then  any  PIO  is  inherently  unacceptable. 

Since  it  is  never  desirable  to  have  an  unstable  closed-loop  pilot/aircraft  system,  the  above 
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definition  dictates  there  can  be  no  acceptable  characteristics  of  a  PIO.  If  there  are  oscillation 
tendencies  due  to  coupling  of  the  aircraft  and  pilot  which  have  relatively  benign  characteristics, 
then  these  should  be  dealt  with  as  separate  handling  qualities  issues.  The  term  PIO  should  be 
reserved  for  description  of  an  unacceptable  coupling  between  the  pilot  and  aircraft  so  that  it  keeps 
its  current  day  connotations  and  continues  to  be  cause  for  the  greatest  concern. 

3.J  Causes 

Several  dynamic  factors  have  been  documented  as  contributing  to  the  occurrence  of  PIO 
These  can  be  separated  into  those  that  are  attributed  to  the  linear  dynamics  of  the  aircraft  and 
pilot,  and  those  that  are  nonlinear  in  nature.  Also,  it  has  been  noted  by  several  researchers  that  a 
trigger  event  is  often  necessary  to  excite  the  PIO.  A  brief  list  of  the  causes  of  PIO  and  a 
description  of  events  associated  as  triggers  is  given  below. 

3.3.1  Linear 

Aspects  of  the  linear  dynamics  of  an  aircraft  that  have  been  associated  with  causing  PIO 
include:  1)  sluggish  response  modes,  2)  lightly  damped  modes,  3)  low  stick  force  per  g,  4) 
sensitive  stick  gradients,  5)  unusual  coupling  responses,  and  6)  excessive  lags  in  the  augmented 
aircraft.  The  most  common  aspect  of  the  linear  pilot  that  is  attributed  to  causing  PIO  is  an 
increased  pilot  gain.  In  the  subsequent  analysis,  the  effect  of  changes  in  pilot  gain  on  system 
stability  is  of  primary  concern 

3.3.2  Nonlinear 

Nonlinearities  associated  with  PIO  include;  l)surface  actuator  rate  and  position  limiting,  2) 
nonlinear  stick  characteristics,  and  3)  mode  switching.  Although  saturation  of  the  control  surface 
actuators  is  common  in  PIOs,  it  is  not  clear  whether  it  is  the  cause  or  effect  of  PIO  The 
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examples  discussed  in  this  report  will  focus  on  the  affect  of  actuator  rate  and  command  limiting  of 
the  augmented  aircraft  on  the  occurrence  of  PIO. 

3.3.3  Trigger  Events 

While  trigger  events  are  not  causes  of  PIO  per  se,  they  are  often  needed  to  excite  a  response 
of  the  pilot  that  will  “trigger”  or  start  the  oscillation  event.  A  trigger  event  can  be  something  as 
simple  as  a  gust  of  wind  or  something  more  complicated  such  as  control  system  mode  switching 
or  a  control  system  failure.  In  the  subsequent  analysis,  closed-loop  stability  is  analyzed  in  the  face 
of  changes  in  both  the  aircraft  and  pilot  models.  That  is,  the  post  trigger  dynamic  situation  is  to 
be  analyzed. 
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4  PIO  Analysis  Process 


The  process  outlined  below  is  intended  to  distinguish  between  aircraft  that  are  PIO  susceptible 
and  aircraft  that  are  PIO  resistant.  The  analysis  is  based  on  determining  the  stability  of  the  closed- 
loop  pilot/aircraft  system  in  the  face  of  variations  in  the  pilot  model  as  well  as  uncertainties  in  key 
aircraft  parameters.  The  process  is  able  to  account  for  many  of  the  nonlinearities  present  in  both 
the  pilot  and  aircraft.  For  the  purposes  of  this  report,  those  analysis  methods  which  include  pilot 
models  are  considered  “closed-loop  methods”  even  if  they  only  consider  the  open-loop 
pilot/vehicle  characteristics.  Methods  which  do  not  include  a  pilot  model  are  considered  “open- 
loop  methods”.  This  point  of  view  is  taken  for  convenience  and  is  not  meant  to  strictly  categorize 
any  particular  method. 

4.1  Process  Outline 

The  underlying  purpose  of  the  analysis  is  to  determine  under  what  conditions  the  pilot/aircraft 
system  is  capable  of  producing  a  sustained  oscillation.  The  unified  multivariable  analysis  method 
consists  of  four  steps.  The  first  step  is  to  determine  an  analytical  model  of  the  human  pilot  that  is 
appropriate  for  the  flying  task  under  study.  Any  accepted  modeling  method  can  be  used,  although 
the  results  may  vary.  Both  crossover  and  "optimal  control"  model  forms  are  acceptable.  A 
synchronous  (constant  gain)  model  is  also  appropriate  for  some  cases  although  a  model  with  at 
least  a  -20  dB/decade  high-frequency  roll-off  characteristic  (to  help  attenuate  feedback  signal 
superharmonics)  is  desirable. 

The  second  analysis  step  is  to  identify  and  isolate  the  dynamic  elements  suspected  of 
contributing  to  PIO  behavior.  Linear  and  nonlinear  elements  can  be  considered  simultaneously. 
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Aerodynamic  nonlinearities,  control  system  nonlinearities,  nonlinear  stick  gradients,  various  pilot 
model  parameters,  and  control  surface  rate  limits  are  of  particular  interest. 

The  third  step  is  to  perform  a  stability  robustness  analysis  using  the  structured  singular  value 
method.  This  analysis  will  yield  a  bound  on  the  largest  allowable  parameter  variation  before 
instability  is  possible.  The  structured  singular  value  analysis  leads  to  a  sufficient  condition  for 
stability  and  may  be  conservative.  If  the  closed-loop  system  operates  in  a  fashion  such  that  the 
predicted  bound  on  the  uncertainty  is  not  exceeded,  then  the  closed-loop  system  will  remain  stable 
and  will  not  limit  cycle.  However,  if  the  predicted  bound  is  exceeded,  the  system  may  or  may  not 
be  unstable  and  may  or  may  not  limit  cycle. 

The  final  analysis  step  is  to  search  for  limit  cycle  behavior  of  the  nonlinear  pilot/aircraft 
system.  A  limit  cycle  is  a  self  sustained  oscillation  that  is  particular  to  nonlinear  systems.  The 
limit  cycle  search  is  carried  out  in  the  frequency  domain.  Results  yield  frequency  and  amplitude  of 
the  oscillation  as  well  as  a  prediction  of  limit  cycle  stability.  Stable  limit  cycles  will  return  to  a 
constant  amplitude  oscillation  following  a  perturbation  whereas  an  unstable  limit  cycle  will  not. 
Limit  cycle  predictions  are  generally  confirmed  using  nonlinear  time  domain  simulations. 

4.2  Open-loop  Analysis 

Though  the  analysis  methods  developed  in  this  report  are  considered  closed-loop,  some 
existing  open-loop  methods  are  mentioned  here  as  they  are  relevant  to  the  construction  of  a 
unified  theory  in  Section  4.3.5.  The  open-loop  characteristic  that  appears  most  closely  related  to 
PIO  susceptibility  is  the  behavior  of  the  aircraft  frequency  response  near  the  pilot/aircraft  gain  and 
phase  crossover  frequencies.  The  Smith-Geddes  Criterion,  for  example,  requires  computation  of 
the  magnitude  slope  of  the  aircraft  frequency  response  in  the  range  between  1 .0  and  6.0 
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rad/sec.  t21]  This  slope  is  then  used  to  determine  a  criterion  frequency  where  the  aircraft  transfer 
function  phase  angle  is  checked.  If  the  phase  angle  is  less  than  -1 80  deg,  a  PIO  is  predicted  and  if 
the  phase  angle  is  less  than  -165  deg,  a  PIO  is  possible.  Nevertheless,  the  criterion  frequency  cocr 
is  defined  as, 

©cr  =  6.0  +  0.248  (4-1) 

where  S  is  the  computed  magnitude  slope  in  dB/octave. 

Since  the  aircraft  transfer  function  phase  angle  becomes  more  negative  at  higher  frequency,  a 
large  negative  amplitude  slope  reduces  the  criterion  frequency  and  naturally  leads  to  a  larger 
phase  angle.  Therefore,  the  Smith-Geddes  PIO  criterion  tends  to  identify  aircraft  configurations 
having  the  characteristic  of  large  negative  magnitude  slopes  in  the  frequency  range  of  1 .0  to  6.0 
rad/s  as  PIO  resistant.  In  addition,  an  aircraft  with  a  shallow  phase  angle  slope  in  the  same 
frequency  range  will  decrease  the  phase  angle  sensitivity  to  criterion  frequency  changes  and  will 
be  predicted  less  PIO  susceptible.  Since  the  Smith-Geddes  Criterion  phase  angle  is  -180  deg,  it  is 
reasonable  to  suggest  that  configurations  with  shallow  phase  angle  slopes  in  this  region  will  be 
predicted  less  PIO  susceptible. 

The  flight  test  data  base  given  in  Reference  [8]  provides  experimental  data  which  can  be  used 
to  evaluate  (longitudinal)  PIO  prediction  methods.  The  data  base  has  been  referred  to  elsewhere 
as  the  HAVE  PIO  data  base,l*  '1  and  that  designation  is  continued  in  this  report.  The  experiments 
were  conducted  using  the  USAF/Calspan  variable  stability  NT-33  aircraft.  The  18  aircraft/flight 
control  system  configurations  were  evaluated  in  landing  approach  tasks.  The  four  baseline 
configurations  (H2-1,  H3-1,  H4-1,  H5-1)  were  chosen  to  have  different  short  period  dynamics 
that  all  met  MIL-F-8785C  Level  1  boundaries  for  Category  C  landing  approach.  The  phugoid 
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and  lateral/directional  characteristics,  which  all  met  Level  1  requirements,  were  held  constant. 
PIO  ratings,  Cooper-Harper  pilot  ratings,  and  pilot  comments  were  obtained. 

For  Configuration  H2-1,  the  control  stick  deflection  to  pitch  attitude  transfer  function  is. 


g(s)  = 


(14) 

s[.64\2.4][.68\26] 


(4-2) 


where  the  shorthand  notation  is  (t)  =  (xs+l)  and  [Ocon]  =  (s/con)^  +  2Cs/(0n  +  1.  The  aircraft 
transfer  function  for  Configuration  H2-5  is. 


g(s)  = 


(1.4) 

s(1.0)[.64\2.4][.68\26] 


(4-3) 


Figure  4-1  shows  a  Bode  diagram  of  g(jco)  for  HAVE  PIO  Configurations  H2-1  and  H2-5.  The 
average  PIO  rating  for  Configuration  H2-1  was  1  0  while  H2-5  was  given  average  rating  of  4.33. 
An  aircraft  receiving  a  rating  above  2.0  is  considered  PIO  susceptible. 

In  the  frequency  range  from  1.0  to  6.0  rad/sec.  Configuration  H2-5  has  a  steeper  magnitude 
slope  than  Configuration  H2-1 ,  but  H2-1  exhibits  a  shallower  phase  angle  slope.  Although  the 
criterion  frequency  for  H2-5  is  lower  than  that  for  H2-1  (approximately  3  rad/sec  as  opposed  to  4 
rad/sec),  the  phase  angle  for  H2-5  is  less  than  -180  deg  and  the  phase  angle  for  H2-1  is  greater 
than  -165  deg  at  the  corresponding  criterion  frequencies.  Therefore,  one  can  see  from  Figure  4-1 
that  Configuration  H2-5  will  be  more  susceptible  to  PIO  than  Configuration  H2-1. 

Mitchell  and  Hoh  advocate  a  PIO  criterion  that  specifies  a  maximum  phase  angle  change  after 


the  phase  angle  crosses  negative  180  deg.1^^1  Their  criterion  is  specified  in  terms  of  a  phase  delay 


parameter,  Xp,  that  is  defined  by. 


A<I> 

^  “  2(57.3) (0,, 


(4-4) 


where  is  the  phase  loss  from  ooigo  to  2Q)igo  and  caigo  is  the  frequency  where  the  aircraft 
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Figure  4-1.  Bode  Plots  for  HAVE  PIO  H2-1  and  H2-5. 

transfer  function  reaches  -180  degrees.  An  aircraft  is  considered  PIO  susceptible  if  the  phase 
delay  Xp  >  .  12  sec  (.  15  sec  in  landing). 

One  can  see  by  (4-4)  that  upper  limits  on  the  phase  delay  parameter  will  tend  to  limit  the 
aircraft  phase  angle  slope  near  the  -180  deg  point.  A  small  value  of  Xp  is  needed  to  categorize  the 
aircraft  as  PIO  resistant,  so  it  will  therefore  require  a  small  value  of  A<I>  and/or  a  large  value  of 
CO  180.  These  requirements,  in  effect,  will  list  aircraft  that  have  shallow  phase  angle  slopes  in  the 
frequency  range  where  the  phase  angle  crosses  -180  deg  as  PIO  resistant. 

In  summary,  these  open-loop  methods  deal  directly  with  the  shape  of  the  aircraft  transfer 
function  near  the  point  where  the  phase  angle  is  -1 80  deg  and  near  where  the  pilot/vehicle 
crossover  frequency  is  likely  to  occur  (1-6  rad/sec). 


15 


4.3  Closed-loop  Analysis 

Tlie  four  steps  in  the  multivariable  PIO  analysis  process  outlined  in  Section  4.1  are  discussed 
in  detail  in  the  following  sections.  Sections  4.3.5  and  4.3.6  use  stability  robustness  ideas  to 
establish  ties  between  existing  methods  and  to  set  forth  a  new  analysis  procedure.  A  detailed 
example  of  the  multivariable  analysis  method  is  then  presented  in  Section  4.3.7. 

4.3.1  Pilot  Modeling 

Advances  in  the  theory  of  pilot  modeling  have  led  to  new  theories  of  pilot-in-the-loop 
oscillations.  Early  PIO  criteria  involved  concepts  from  the  crossover  model  of  the  human 
operator  while  more  recent  research  has  focused  on  the  optimal  control  model.  Kleinman's 
optimal  control  model  (OCM)  of  the  human  operator  is  based  upon  Linear,  Quadratic,  Gaussian 
(LQG)  control  theory.  The  model  includes  a  set  of  state-feedback  gains  that  are  computed  using 
a  special  quadratic  cost  function.  The  state-feedback  control  law  is  then  augmented  with  a  state 
estimator  and  a  predictor.  The  purpose  of  the  predictor  is  to  model  the  human's  ability  to  adapt  to 
his/her  inherent  time  delay. 

A  modified  version  of  the  Kleinman  model  has  been  developed  by  Davidson  and  Schmidt. l^-^l 
The  modified  optimal  control  model  (MOCM)  replaces  the  pure  time  delay  with  a  Pade' 
approximation  of  the  time  delay.  The  Pade'  approximation  used  by  Davidson  and  Schmidt  is 
given  by, 

(4-5) 

l  +  r,,,v/2+(r,,.v)  /x 

where  Tj  is  the  human  pilot  time  delay.  The  principal  advantage  of  the  Pade'  approximation  is  that 
its  use  allows  for  direct  computation  of  a  human  operator  transfer  function.  Whereas  Kleinman's 
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model  yields  only  frequency  response  (describing  function)  information  of  the  human  response, 
the  Davidson  model  provides  detailed  information  of  dominant  response  modes.  Therefore,  the 
interaction  between  the  pilot  and  aircraft  can  be  more  readily  identified. 

Anderson  has  recently  revised  the  Davidson  model  to  use  state-space  solutions  for  both  H2 
and  Hjkj  norm  minimizing  controllers. l^l  This  latest  extension  uses  the  "standard"  optimal  control 
solutions  that  have  been  published  recently  by  Doyle,  et.  al.l^^l  The  "standard"  model  form  is 
given  by  the  following  state-space  model. 


For  application  to  pilot  modeling,  the  vector  Uj  represents  the  external  noise  or  command,  U2  is 
the  vector  containing  the  pilot  manipulator  input,  yj  is  the  vector  representing  the  pilot's  task 
objective,  and  represents  the  pilot's  observation  variables. 

Because  the  newer  “standard”  H2/H„  model  has  not  yet  been  validated  for  aircraft 

applications,  Davidson  and  Schmidt’s  MOCM  is  used  throughout  this  work.  The  specific  MOCM 
pilot  models  used  in  the  examples  presented  in  this  report  share  the  modeling  parameters  shown  in 
Table  4-1. 

Table  4-1.  MOCM  Pilot  Model  Parameters. 


Parameter _ _ _ Value 

time  delay  0  2  sec 

neuromuscular  time  constant  0. 1  sec 

observation  signal-to-noise  ratio  0.01 

neuromotor  signal-to-noise  ratio  0.003 
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The  optimal  MOCM  cost  function  is  assumed  to  be  given  by, 

J  =  ^jp(t)  +  P5:{t)]dt  (4-7) 

where  e(t)  is  the  error  signal  and  6y(t)  is  the  pilot  manipulator  command  input.  The  control  rate 

weighting  p  is  chosen  to  achieve  the  desired  neuromuscular  time  constant  as  described  in 
Reference  [14],  Note  that  the  error  rate  is  retained  in  the  observation  vector  but  is  not  weighted 
in  the  cost  function  (4-7). 

Smith  has  described  a  contactor  pilot  model  for  study  of  fully  developed  PIO  behavior.  1^51 
The  contactor  nonlinearity  is  essentially  a  relay  that  yields  a  fixed  positive  value  for  any  positive 
value  input  and  a  fixed  negative  value  of  any  negative  value  input.  For  example,  if  the  input  to  the 
nonlinearity  is  +1,  the  relay  output  might  be  +2.  For  an  input  of  +4,  the  output  of  the  same  relay 
is  still  +2.  This  model  will  also  be  used  in  Section  4. 3. 5. 4. 

4.3.2  Identify  Dynamic  Elements  Suspected  of  Contributing  to  PIO 

With  given  models  of  the  aircraft  dynamics  and  the  human  pilot,  the  elements  suspected  of 
contributing  to  PIO  must  be  identified  and  isolated.  The  elements  directly  associated  with  the 
causes  of  PIO  given  in  Section  3.3  can  be  considered  here.  Since  the  pilot  model  plays  a  key  role 
in  the  current  analysis,  any  of  the  pilot  model  parameters  may  need  to  be  considered  as  uncertain 
elements  that  can  contribute  to  PIO  tendency. 

4.3.3  Stability  Robustness  Analysis 

A  PIO  is  fundamentally  a  stability  problem.  Sustained  oscillations  only  occur  in  linear  systems 
when  the  dynamics  of  the  system  lose  damping.  An  example  is  an  undamped  spring-mass  system 
that  will  oscillate  from  an  initial  displacement  with  constant  amplitude  motion. 
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Stability  robustness  analysis  is  the  study  of  how  model  parameter  variation  or  "uncertainty" 
affects  feedback  loop  stability.  Since  the  mid  1980's,  a  resurgence  in  frequency-domain  analysis 
methods  for  stability  robustness  has  occurred.  These  methods  have  relied  on  the  singular  value 
matrix  norm  as  a  "measure"  of  relative  stability.  For  example,  gain  and  phase  margins  for 
multivariable  feedback  systems  can  be  defined  using  the  singular-value  based  analysis  method.  I2<>l 
Stability  margins  defined  in  such  a  manner  are  very  useful  in  that  they  define  the  gain  or  phase 
variation  that  can  occur  simultaneously  in  each  feedback  loop.  As  a  result,  an  investigation  of 
multiple  parameter  variations  can  be  carried  out  rather  easily. 

Doyle  is  widely  regarded  as  the  key  developer  of  the  structured  singular  value  analysis 
method.P'^]  The  structured  singular  value  method  allows  the  research  engineer  to  investigate  not 
only  simultaneous  model  parameter  variations,  but  also  the  effect  of  multiple  parameters  located 
arbitrarily  throughout  the  feedback  system.  Thus,  parameter  variation  is  not  confined  to  just  the 
feedback  loop  signals,  but  can  include  individual  system  parameters. 

4,3.3. 1  Basic  Theory 

Perhaps  the  most  important  contribution  of  the  structured  singular  value  method  is  the 
construction  of  analysis  block  diagrams  like  Figure  4-2.  The  transfer  function  matrix  M(s)  in 
Figure  4-2  does  not  necessarily  represent  the  aircraft,  the  pilot,  or  a  series  combination  of 
feedback  loop  elements.  The  function  M(s)  represents  the  transfer  function  matrix  as  seen  by 
each  model  parameter  that  is  varying  -  regardless  of  where  the  parameter  occurs  in  the  feedback 
loop.  For  example,  a  variation  due  to  inaccurate  aerodynamic  modeling  can  be  analyzed  in 
conjunction  with  a  variation  in  the  pilot's  loop  gain  Consequently,  model  variations  from  many 
different  sources  can  be  analyzed  at  once. 
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Figure  4-2.  Stability  Robustness  Analysis  Diagram. 

Stability  of  the  feedback  loop  can  be  investigated  from  the  loop  equations  represented  in 
Figure  4-2, 

y(s)  =  -A  x(s)  (4-8) 

x(s)  =  M(s)  y(s)  (4-9) 

Solving  (4-9)  in  terms  of  x(s)  yields, 

[I  +  M(s)A]  x(s)  =  0  (4-10) 

The  matrix  [I  +  M(s)A]  is  known  as  the  return  difference  matrix.  Stability  robustness  analysis 
involves  finding  the  size  of  the  (block)  diagonal  matrix  A  that  satisfies  det[I  +  M(ja))A]  =  0. 

Doyle,  therefore,  defined  the  structured  singular  value  |J.(M,(o)  as, 

— =T  {amax(A)|det[I  +  M0a))A]  =  O}  (4-11) 

where  Omax  ‘s  the  largest  matrix  singular  value.  The  structured  singular  value  computation 
attempts  to  find  the  smallest  uncertainty  block  A  that  can  destabilize  the  feedback  system.  Since 
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the  minimization  problem  of  Equation  (4-11)  can  be  difficult  to  solve,  Doyle  derived  expressions 
for  upper  and  lower  bounds  on  the  structured  singular  value.  Namely, 

p(UM)  <  f.i(M,a))  <  amax(DMD‘')  (4-12) 

U  D 

where  U  is  unitary  and  D  is  a  real  scaling  matrix  such  that  DAD''  =  A.  The  lower  bound  is  not 
used  as  a  measure  of  robustness  but  is  used  to  determine  how  close  the  upper  bound 
approximates  the  exact  value.  For  practical  systems  the  upper  and  lower  bounds  differ  by  less 
than  15%  showing  the  upper  bound  to  be  a  reasonable  approximation.l^S]  The  upper  bound  is  still 
posed  as  an  optimization  problem,  but  it  is  one  that  has  no  non-global  minima.  This 
optimization  involves  finding  the  scaling  matrix  that  minimizes  the  maximum  singular  value  of 
DMD  '.  The  optimum  scaling  matrix  D  can  be  found  through  standard  nonlinear  optimization 
techniques,  but  is  often  approximated  using  an  iterative  method  first  developed  by  Osborne.l^*4 
Osborne’s  technique  finds  that  D  which  minimizes  the  Frobenius  norm  of  DMD  ‘.  Because  the 
Frobenius  norm  acts  as  an  upper  bound  on  the  2-norm  (maximum  singular  value),  reducing  the 
Frobenius  norm  helps  reduce  the  2-norm.  P^l  Since  the  Frobenius  norm  is  an  upper  bound  on  the 
2-norm,  use  of  Osborne’s  technique  to  select  D  yields  a  conservative  estimate  of  the  stability 
margin,  but  the  technique  is  very  efficient. 

As  given  above,  the  elements  of  A  can  be  complex.  When  the  parameter  variations  modeled  in 
A  are  known  to  be  real  only,  a  less  conservative  upper  bound,  developed  by  Jones  128],  can  be 
found  for  the  structured  singular  value.  This  upper  bound  is  given  by, 

t4(M,(o)<'^"  amax[(DM0co)D-la>  +  (DM0co)D-l<D)T)/2]}  (4-13) 

where  D  is  determined  as  before.  The  set  of  matrices  defined  by  <!>  contains  all  possible 
combinations  of +1  on  the  diagonal.  These  matrices  represent  the  possible  directions  of  variation 
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(positive  or  negative)  of  each  of  the  real  uncertainties.  Use  of  the  upper  bound  for  real  parameter 
variations  yields  not  only  the  minimum  size  of  A  that  can  destabilize  the  system,  but  also  the  worst 
case  directions  for  each  of  the  individual  variations. 

The  structured  singular  value  analysis  method  can  also  be  used  to  study  the  effect  of  isolated, 
single-valued  nonlinear  elements  on  system  stability.  If  a  sector  can  be  drawn  around  the  graph  of 
a  nonlinear  element,  then  the  structured  singular  value  can  be  used  to  determine  if  the  sector- 
bounded  nonlinearity  can  destabilize  the  closed-loop  system.  References  [31]  and  [32]  describe 
analysis  of  nonlinear  saturation  using  structured  singular  value  analysis  and  sector-bounding 
methods,  respectively. 

Whether  the  matrix  A  represents  linear  or  sector-bounded  nonlinear  elements,  the  result  of  the 
structured  singular  value  analysis  is  only  the  minimum  size  of  A.  The  character  of  the  instability 
that  may  result  if  the  size  of  A  is  greater  than  the  minimum  is  not  known.  For  example,  it  is  not 
known  whether  the  resulting  instability  is  purely  exponential  or  oscillatory. 

4. 3. 3, 2  M-A  Construction 

An  essential  part  of  the  stability  robustness  analysis  is  the  construction  of  the  M  matrix  shown 
in  Figure  4-2.  A  first  step  involves  identification  of  the  uncertainties  to  be  analyzed.  Three  basic 
types  of  uncertainty  will  be  considered  herein.  Figure  4-3  shows  a  multiplicative  uncertainty  and  a 
additive  uncertainty.  The  variation  (5)  represents  a  percentage  signal  change  for  the  multiplicative 
uncertainty  and  an  absolute  change  for  the  additive  uncertainty.  Figure  4-4  illustrates  the 
modeling  of  a  nonlinear  saturation  by  an  equivalent,  but  uncertain  gain.  The  effect  of  a  saturation 
element  is  to  reduce  the  effective  gain,  down  from  unity,  when  the  input  to  the  nonlinear  element 
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a)  Multiplicative  Uncertainty 


b)  Additive  Uncertainty 

Figure  4-3,  Representation  of  Uncertainty. 

Saturation  Equivalent  Gain  Representation 

Figure  4-4.  Representation  of  a  Nonlinear  Saturation  Element  by  an  Equivalent  but 
Uncertain  Gain. 

exceeds  the  saturation  limit.  Therefore,  the  effective  gain  is  a  fiinction  of  the  input  to  the 
nonlinear  element. 

Once  the  uncertainties  have  been  identified,  the  system  must  be  put  in  the  standard  M-A  form 
for  analysis.  While  some  algorithms  exist  to  aid  in  the  construction  of  the  necessary  matrix,  l  ^^l 
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it  is  suggested  that  a  graphical  simulation  program  like  SIMULINK’^”  be  used  to  determine  the 
state-space  model  for  M(s).  SIMULINK™  allows  the  user  to  linearize  a  feedback  system  and 
define  the  input  and  output  from  any  signal  in  the  feedback  loop.  To  linearize  a  system, 
SIMULINK^^  inport  and  outport  elements  are  placed  around  the  uncertainties  as  illustrated  in 
Figure  4-5.  This  diagram  is  for  multiplicative  uncertainty  shown  in  Figure  4-3(a).  Note  a 
negative  sum  is  needed  to  generate  the  proper  representation  of  M(s)  as  defined  in  Figure  4-2. 
With  all  the  inport/outport  pairs  in  place,  a  state-space  model  of  M(s)  is  produced  using  the 
MATLAB^“  linmod  command.  With  M(s)  computed,  the  affect  of  the  modeled  uncertainties  on 
system  stability  is  analyzed  using  the  structured  singular  value  method.  Note  that  the  system 
under  study  can  be  completely  linear  (and  does  not  need  to  be  linearized).  However,  the  linmod 
feature  of  MATLAB^'^  provides  a  convenient  and  useful  way  to  form  M(s)  for  complicated  block 
diagrams. 


Figure  4-5.  Model  of  Uncertainty  in  a  SIMULINK™  Diagram. 

4.3.3.3  Numerical  Algorithms 

MATLAB™  contains  two  routines  (in  the  Robust  Control  Systems  toolbox)  to  calculate  the 
upper  bound  on  the  structured  singular  value  {ssvbode  and  perron).  Both  of  these  algorithms 
compute  the  structured  singular  value  bound  of  Equation  (4-12)  which  considers  complex 
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uncertainties.  These  algorithms  are  based  on  the  work  of  Safonov,  These  routines  can  yield 

overly  conservative  estimates  for  systems  with  only  real  parameter  variations. 

For  the  immediate  purpose  of  analyzing  pilot/aircraft  systems,  only  real  uncertainties  have 

been  considered.  Therefore,  an  algorithm  based  on  that  outlined  by  Reference  [29]  has  been 

written  to  calculate  the  bound  of  Equation  (4-13).  The  algorithm  outline  is  given  below. 

Algorithm  to  Compute  Upper  Bound  on  the 
StnirtiireH  Singular  Value  for  Real  Parameter  Variations 

1 .  Define  frequency  range  of  interest,  (Oi  <  ca  <  (Of. 

2.  Set  up  2‘’x  q  matrix  (Oj)  whose  rows  account  for  the  possible  permutations 
of  1  and  -1  (q  =  number  of  uncertainties). 

3.  Over  all  co, 

a)  Calculate  the  transfer  function  matrix  at  s  =  jeo,  i.c.  M(iw). 

b)  Use  Osborne’s  technique  to  find  D. 

c)  For  each  row ,  k ,  in  Ot  , 

1)  Calculate  diagonal  0  matrix  as  <b  =  diag[k‘'*  row  of  (DtI 

2)  Find  a™ax[(DM(jco)D-‘0  +  (DM(j«))D-ld))T)/2]. 

Store  result  as  SV(k). 

d)  Find  the  maximum  value  in  SV(k)  calculated  above. 

This  represents  |Jt(M,(o)  at  this  frequency. 

4.  Plot  the  log(pr(M,(o))  versus  log(co). 

The  maximum  (real)  structured  singular  value  over  all  frequencies  will  be  denoted  Pr  and  will 
be  used  to  define  the  maximum  allowable  variation  while  guaranteed  stability  is  maintained. 
Assuming  D=diag[6i,62,...,6q],  the  closed-loop  system  will  remain  stable  as  long  as  |6i|  <  1/pr. 

4.3.3.4  Verification 

To  ensure  the  stability  robustness  algorithms  work  as  desired,  an  analysis  of  two 
configurations  of  the  HAVE  PIO  database  was  conducted.  Configurations  H2-1  and  H2-5  were 
modified  to  include  a  first  order  actuator  as  shown  in  Figure  4-6.  Table  4-2  gives  the  vehicle  and 
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Figure  4-6.  HAVE  PIO  Configuration  Model  Modified  to  Include  Actuator  Dynamics, 


Table  4-2.  Vehicle  and  Pilot  Model  Parameters  for  Modified  HAVE  PIO  Configurations. 


®sp  / 

T2 

Xr 

Tpi 

Tp2 

T 

Configuration 

(rad/s/-) 

(s) 

(s) 

L 

Kp 

(S) 

(S) 

(s) 

H2-1 

2.4/0,64 

0 

0.05 

1.0 

0.74 

0.27 

0 

0.3 

H2-5 

2.4/0.64 

1.0 

0.05 

1.0 

0.52 

2.37 

0 

0.3 

pilot  model  parameters  used  in  the  analysis.  The  pilot  model  parameters  were  chosen  using  rules 
from  the  Neal-Smith  criteria,  l^^l 

The  two  dynamic  effects  to  be  considered  here  are  a  simultaneous  change  in  pilot  gain  and 
control  surface  actuator  rate  limiting.  The  rate  limit  is  modeled  by  a  saturation  element  in  the 
actuator  model  and  is  represented  by  an  equivalent  but  uncertain  gain.  A  representation  of  the 
SIMULINK’^^  diagram  used  to  define  M(s)  is  shown  in  Figure  4-7.  The  first  variation  represents 
a  percentage  change  in  the  pilot  gain  while  the  second  represents  a  percentage  change  in  the  unity 
gain  model  of  the  saturation  element. 

Figure  4-8  shows  the  inverse  real  structured  singular  value  calculation  for  both  configurations 
using  Equation  (4-13).  The  minimum  of  the  inverse  structured  singular  value  plot  reveals  the 
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M 


Figure  4-7.  SIMULINK™  Diagram  Used  to  Define  M(s)  for  Structured  Singular  Value  Analysis. 


Frequency  (rad/s) 


Figure  4-8.  Inverse  Structured  Singular  Value  of  M(s). 

largest  bound  on  5]  and  62  (from  Figure  4-7  )  before  instability  is  possible.  The  minimum  value  is 
0.56  for  configuration  H2-1  and  0.08  for  configuration  H2-5.  The  structured  singular  value 
calculation  of  Equation  (4-13)  also  reveals  the  direction  of  the  worst  case  variation.  At  the 
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frequency  where  the  minimum  is  found,  the  matrix  <I)  that  maximized  (4-13)  is  <I)  =  diag[+l,  -1], 
Therefore,  the  worst  case  variation  in  pilot  gain  is  an  increase  of  (1+  6|)  and  the  worst  case  gain 

variation  due  to  rate  limiting  is  (1-  52). 

The  stability  boundary  calculated  from  the  real  structured  singular  value  analysis  of 
configurations  H2-1  and  H2-5  are  shown  in  Figures  4-9  and  4-10  respectively.  The  shaded  blocks 
represent  values  for  pilot  gain  and  equivalent  gain  due  to  rate  limiting  for  which  the  structured 
singular  value  analysis  guarantees  stability.  A  direct  calculation  of  the  system  poles  yields  the 
stability  boundaries  shown  by  the  solid  curves.  The  dashed  lines  represent  the  nominal  values  of 
the  pilot  gains  in  each  case.  The  nominal  value  for  the  equivalent  gain  due  to  rate  limiting  is  L  = 

1.  Figures  4-9  and  4-10  illustrate  how  well  the  structured  singular  value  approximation  fills  the 
stable  region  in  the  parameter  space.  The  corners  of  the  structured  singular  value  boundaries 
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Figure  4-9.  Stability  Boundary  For  Configuration  H2-1 . 
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Figure  4-10.  Stability  Boundary  for  Configuration  H2-5. 

closely  approach  the  curves  separating  the  stable  and  unstable  regions  in  the  parameter  space 
indicating  the  structured  singular  value  to  be  a  good  approximation  to  the  actual  stability 
boundary. 

4.3.3.S  Summary 

When  A=0,  M(s)  represents  the  nominal  system.  When  the  uncertain  parameters  are  known  to 
vary  over  a  given  range,  they  are  usually  scaled  such  that  -1  <  5  <  1  encompasses  the  given 
variation.  With  a  scaled  system  and  l/pr  ^  1,  the  system  can  handle  the  given  variation  and 
remain  stable.  If  1/pr  ^  1,  then  the  system  can  tolerate  variations  up  to  at  least  l/pr  while 
maintaining  system  stability.  For  the  case  of  a  variation  in  a  unity  gain  representing  a  saturation 
element,  the  value  of  l/pr  represents  the  percent  saturation  of  the  element.  That  is,  the  system  can 
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tolerate  an  input  to  a  nonlinear  element  up  to  Pr  /(1-Pr)  percent  greater  than  the  saturation  limit 
while  guaranteed  stability  is  maintained. 

4.3.4  Multivariable  Describing  Functions 

Some  researchers  have  noted  that  PIOs  are  analogous  to  limit  cycles.l'^^^’^^l  A  limit  cycle  is  an 
unforced,  sustained  oscillation  of  a  nonlinear  system.  The  dynamics  associated  with  pilot/aircraft 
interaction  are  fundamentally  nonlinear.  In  fact,  both  the  pilot  and  aircraft  are  very  complicated 
nonlinear  systems.  Graham  and  McRuer  also  note  that  friction  and  hysteresis  characteristics  of 
early  jet  fighter  flight  control  systems  led  to  serious  PlOs.t^^’l  Thus,  methods  for  predicting  limit 
cycles  for  nonlinear  systems  should  also  be  considered  in  the  development  of  new,  unified  PIO 
analysis  methods. 

A  useful  class  of  nonlinear  systems  includes  nonlinearities  isolated  as  separate  entities  while 
the  remaining  linear  elements  are  represented  in  lumped-parameter  form.  For  this  particular  class 
of  systems,  an  analysis  based  upon  quasi-linear  harmonic  approximations  is  readily  completed. 

The  quasi-linear  approximations  are  usually  known  as  describing  functions,  For  PIO 
analysis,  a  describing  function  approach  is  particularly  attractive  as  it  has  close  ties  with  linear 
frequency-domain  analysis.  Most  existing  PIO  analysis  methods  are  also  based  on  control- 
theoretic,  frequency-domain  methods. 

4.3.4. 1  Basic  Theory 

Assume  the  input  to  a  single-input,  single-output  nonlinear  element  is  sinusoidal  with  a  given 
amplitude  and  frequency.  The  output  of  the  nonlinear  element  will  generally  not  be  sinusoidal, 
but  will  be  periodic  with  the  same  period  as  the  input  sinusoid.  If  the  output  wave  is  then 
analyzed  in  terms  of  its  Fourier  components,  the  fundamental  component  can  be  related  to  the 
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input  in  familiar  terms  of  an  amplitude  ratio  and  a  phase  angle.  The  complex  ratio  of  the 
fundamental  component  of  the  output  to  the  input  is  defined  as  the  sinusoidal  input  describing 
function  (SIDF).!^’^] 

For  a  sinusoidal  input  to  the  nonlinear  element,  x(t)  =  a  sin((»t),  the  output  y(t)  is  expressed  in 
terms  of  its  Fourier  series.  Namely, 

no 

y(t)  =  aj  +  X  (a„cos  ncot  +  b„sin  not)  (4-14) 

n»l 

where, 

a„ -if  y(t)cosno)td((Bt) 

=  if  y(t)sinnratd(a)t) 

For  a  skew  symmetric  nonlinearity,  ao=  0,  the  fimdamental  component  of  the  output  becomes, 

y,(t)  =  a,  coswt  +  b,  sincot  (4-15) 

The  describing  function  is  then  given  by, 

N(a,{D)  =  ^Z(|),  ,  <j),  =  tan''^-^j  (4-16) 

and  a|  and  bj  are  obtained  using  the  integral  expressions  in  (4-14).  In  general  the  describing 
function  can  depend  on  both  the  amplitude  and  frequency  of  the  input  signal.  However,  when  the 
nonlinearity  contains  no  energy  storage  element,  N  is  a  function  of  amplitude  only.t-^^1  All  simple 
(single  valued)  nonlinearities  and  many  complex  nonlinearities  have  describing  fiinctions  which  are 
frequency  invariant.  Table  4-3  lists  4  nonlinearities  and  their  frequency  invariant  describing 
functions. 
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Table  4-3.  Four  Nonlinearities  and  the  Corresponding  Sinusoidal  Input  Describing  Functions. 


Nonlinear  Element 

Sinusoidal  Input  E)cscribing  Function 

gB| 

2m 

7t 

m 

(for  a/a  ^1) 

B 

sin  ^(a/a)  +  (a/a)yi-(a/a) 

(for  a/a  >1) 

n 

0 

(for  a/a  <1) 

-a 

2ni 

sin~^(a/a)  +  (a/a)^l  -  (a/a)^ 

(for  a/a  >1) 

/ 

a 

m - 

7C 

4a 

■ 

- 7 

*a 

Tia 

4a 

m  -1 _ 

B 

■ 

111  1 

Tia 

Figure  4-11  shows  an  example  of  replacing  a  nonlinear  element  by  its  SIDF  resulting  in  a 
quasi-linear  feedback  system.  The  remnant  consists  of  the  higher  harmonics  of  the  Fourier  series 
of  the  output  (terms  with  n  >  1  in  Equation  (4-14)).  The  output  of  the  actual  nonlinearity  is  equal 
to  the  sum  of  the  input  times  the  SIDF  plus  the  remnant.  In  considering  only  the  fundamental 
harmonic  when  developing  the  describing  function,  it  is  automatically  assumed  that  the  remnant 
has  little  affect  on  the  system  output.  This  assumption  holds  if  the  linear  part  of  the  system  is 
sufficiently  low-pass  as  to  attenuate  the  higher  harmonics  present  in  the  remnant.  Also,  the 
maximum  amplitude  present  in  the  remnant  is  normally  much  less  than  the  magnitude  of  the 
fundamental  harmonic.  For  the  case  of  a  saturation  nonlinearity  and  a  K/s  type  plant.  Reference 
[37]  states  the  largest  term  in  the  remnant  (which  is  the  third  harmonic)  is  never  more  than  28% 
of  the  amplitude  of  the  fundamental  component  of  the  output  of  the  nonlinearity.  When 
attenuated  by  the  integrator,  the  amount  of  the  third  harmonic  present  in  the  error  signal ,  x(t),  as 


32 


a)  Nonlinear  System 


Remnant 


b)  Quasi-Linear  System 

Figure  4-11.  Replacement  of  a  Nonlinear  Element  by  a  Sinusoidal  Input  Describing  Fimction. 

compared  to  the  fundamental  is  not  more  than  9%.  This  value  is  a  measure  of  the  error 
introduced  by  ignoring  the  contribution  of  the  remnant. 

Although  the  above  discussion  considers  a  single-loop  system  with  a  single  nonlinearity, 
extension  to  a  multi-loop  multiple-nonlinearity  system  is  straightforward.  Let  A  from  Figure  4-2 
represent  a  diagonal  matrix  of  nonlinear  elements  and  M(s)  now  represents  the  linear  part  of  the 
system  as  seen  by  the  nonlinear  elements  in  A.  Assume  that  the  input  signal  to  each  nonlinear 
element  in  A  is  sinusoidal  and  oscillating  at  the  same  frequency  ca,  so  that, 

xi  =  ai  sin(a»t  +  0i)  (4-17) 

Again,  to  use  describing  function  approximations,  the  linear  part  of  the  system  represented  by 
M(s)  is  assumed  to  act  like  a  low-pass  filter  so  that  superharmonic  signal  components  can  be 


neglected  in  the  analysis.  Under  these  assumptions,  the  nonlinear  elements  can  be  replaced  by 
their  sinusoidal-input  describing  functions  (SIDF),  i.e., 

A^N(A,a))  (4-18) 

where  'A'  is  a  vector  containing  the  amplitudes  of  the  sinusoids  entering  each  nonlinear  element 
and  N(A,a))  is  a  diagonal  matrix  of  describing  functions  made  up  of  individual  elements  like  those 
in  Table  4-3. 

In  order  to  investigate  the  stability  of  the  multivariable  nonlinear  feedback  system,  solution  of 
the  system  characteristic  equation  is  considered.  Since  it  is  assumed  that  each  feedback  loop  is 
oscillating  at  the  same  frequency,  the  Laplace  operator  in  (4-10)  can  be  replaced  with  s=jco. 
Equation  (4- 1 0)  then  becomes, 

[I  +  M(j(D)N(A,a))]  x(jto)  =  0  (4-19) 

Equation  (4-19)  represents  the  harmonic-balance  equations  for  the  system.  Solutions  to  (4-19) 
predict  the  existence  of  a  limit  cycle. 

4, 3. 4. 2  M-A  Construction 

Within  the  current  framework,  it  is  convenient  to  cast  the  multivariable  describing  function 
problem  in  a  form  similar  to  the  stability  robustness  analysis.  Using  a  standard  M-A  configuration 
to  define  the  structure  of  the  nonlinear  system  brings  a  needed  organization  to  the  general 
problem  of  predicting  the  occurrence  of  limit  cycle  behavior.  In  this  case,  M(s)  represents  the 
transfer  function  matrix  as  seen  by  each  nonlinear  element  and  as  such  can  contain  parts  of  the 
vehicle,  control  system,  and  pilot.  It  is  suggested  that  a  state-space  representation  of  M(s)  be 
found  in  a  manner  analogous  to  that  of  Section  4. 3. 3.2.  Figure  4-12  shows  how  the 
SIMULINK^“  inport  and  outport  blocks  are  used  to  separate  out  a  typical  nonlinear  element. 
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Note  a  gain  of-1  must  be  added  at  either  the  inport  or  outport  if  the  form  of  Equation  (4-19)  is  to 
be  used.  If  no  such  gain  is  added,  then  the  negative  of  the  transfer  function  matrix  (-M(s))  is 
returned  from  the  linmod  command. 


1 


Nonlinearity 


Outport 


Inport 


■> 


Figure  4-12.  Model  of  Nonlinear  Element  in  SIMULINK™  Diagram  for  M-A  Construction. 


4.3.4.3  Numerical  Algorithms 

Since  x(j(D)  =  0  represents  a  trivial  solution  to  (4-19),  many  researchers  have  used  the 
following  determinant  to  predict  limit  cycles, 

det  [I  +  M0(o)N(A,(d)]  =  0  (4-20) 

Alternatively,  one  could  search  for  values  of  co  and  A  that  yield  one  or  more  zero  eigenvalues  of 
[I  +  M(jco)N(A,©)]  or  eigenvalues  of  [M(i(o)N(A,oo)]  equal  to  -l+jO.l'*'^! 

There  are  several  computational  methods  available  to  aid  in  the  search  for  limit  cycle  behavior. 
Besides  searching  (4-19)  or  (4-20)  directly,  Gray  and  Nakhia  use  Gershgorin  bands  as  an 
eigenvalue  check.  I'*  Gray  and  Taylor  have  developed  a  search  method  that  considers  each 
feedback  loop  in  a  sequential  order,  while  systematically  eliminating  solution  possibilities 
throughout  the  search.l'’^!  Chang,  et.  al.  have  recently  extended  the  sequential  loop  method  of 
Gray  and  Taylor  to  include  a  homotopy  algorithm  with  claims  of  global  convergence.  ^t 
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about  the  same  time,  Pillai  and  Nelson  reported  search  method  results  based  upon  Newton- 
Raphson  optimization  methods.H^l 

The  method  of  Chang,  et.  al.  warrants  further  explanation  as  it  is  easily  implemented  under  the 
current  M-A  formulation.  Also,  the  use  of  a  linear  homotopy  to  improve  the  convergence  of  the 
search  algorithm  is  desirable.  When  cast  in  the  standard  M-A  form  above,  the  set  of  nonlinear 
equations  that  Chang  refers  to  as  “auxiliary  characteristic  equations”  reduce  exactly  to  equation 
(4-19).  This  equation  can  be  rewritten  as  a  set  of  n  complex  nonlinear  equations  (one  for  each 
nonlinear  element)  of  the  form, 

fj  =  x.(j(o)  +  jMi^(ja))N,(A,,a))x,(ja))  =  0  (4-21) 

k=l 

where  Xi(ja))  is  the  result  of  replacing  s=ja)  in  the  Laplace  transform  of  Xi(t)  of  Equation  (4-17). 
With  some  straightforward  substitutions  and  cancellations,  it  can  be  shown  that  for  numerical 
solution,  Xifjo))  in  the  above  equation  can  be  replaced  by, 

XjOffl) ->  x(Aj,0i)  =  Aje^®-  (4-22) 

The  existence  of  a  limit  cycle  is  then  found  by  solving  the  system  of  n  complex  equations  given  in 
(4-21),  with  the  substitution  of  (4-22),  for  the  2n  unknowns  (Ai,a),A2,02,.  .  ,A„,6n)  The  phase 
angle  of  the  sinusoid  entering  the  first  nonlinear  element  is  taken  to  be  zero  (0i=O)  without  loss  of 
generality. 

When  a  good  initial  guess  of  the  solution  is  available,  the  system  of  nonlinear  equations  can 
likely  be  solved  by  any  method  normally  used  for  nonlinear  systems  (e  g.  Newton-Raphson). 
However,  since  an  adequately  good  guess  is  not  always  available,  a  homotopy  method  of  solution 
is  recommended.  Define  the  linear  homotopy  H  as, 

H(z,  t)  =  ( 1  -  t)(z  -  z„ )  + 1  f(z)  (4-23 ) 
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The  solution  of  H  =  0  for  allt  e  [0,l]  forms  a  space  trajectory.  For  t  =  0,  the  solution  is  simply 


z  =  z„ .  At  t  =  1,  the  solution  is  that  desired  H  =  f  =  0 .  Therefore  if  one  can  trace  the  solution 
H  =  0  from  t  =  0  to  t  =1,  the  desired  solution  will  have  been  reached.  One  way  to  trace  the 
solution  is  to  integrate  the  differential  equation  obtained  by  taking  the  total  derivative  of  H  with 
respect  to  t.  Namely, 

dz 


aH(z,t) 

dz 


aH(z,t) 

dt 


(4-24) 


where, 


'dr, 

dr, 

0H(z,t) 

=  (i-t)i  +  t^^ 

cf(z)_ 

dz, 

dz 

dz 

dz 

i  . 

dz. 

aH(z,t) 

dt 

=  -(z  -  z„)  +  f(z) 

The  necessary  partial  derivatives  can  be  computed  numerically  if  desired. 

When  a  solution  is  obtained,  it  indicates  the  frequency  of  the  resulting  oscillation  along  with 
the  input  amplitudes  and  phase  angles  to  each  isolated  nonlinear  element.  Since  the  solution  of  a 
limit  cycle  governs  the  behavior  of  a  dynamic  system,  it  is  reasonable  to  speak  of  the  stability  of 
the  limit  cycle.  When  a  system  that  has  fallen  into  a  stable  limit  cycle  is  perturbed  slightly  from 
the  limit  cycle  condition,  the  system  will  return  to  the  limit  cycle.  If  a  system  were  able  to  operate 
at  an  unstable  limit  cycle  condition  and  was  perturbed  slightly  from  that  condition,  it  would  not 
return.  The  resulting  behavior  could  be  divergent,  asymptotically  stable,  or  could  converge  to  a 
different  (stable)  limit  cycle.  For  this  reason,  an  unstable  limit  cycle  cannot  be  observed 
experimentally.  According  to  Changl^^]^  the  stability  of  the  limit  cycle  is  determined  by 
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considering  a  real  perturbation  in  the  limit  cycle  frequency,  s  =  a  +  j(o.  The  limit  cycle  will  be 
stable  as  long  as  da/dAi  <  0  . 

Grahaml-^'^l  further  distinguishes  between  soft  and  hard  self  excitation  of  a  stable  limit  cycle. 
Soft  self  excitation  exists  when  any  nonzero  initial  condition  will  result  in  the  sustained  oscillation 
of  the  stable  limit  cycle.  Soft  self  excitation  is  common  in  systems  where  the  linear  transfer 
function  matrix  describes  an  unstable  system.  On  the  other  hand,  hard  self  excitation  exists  when 
the  input  and  subsequent  removal  of  a  certain  amplitude  or  wave  form  signal  is  required  to  excite 
the  oscillation. 

4. 3. 4. 4  Verification 

To  ensure  the  multivariable  describing  function  analysis  methods  work  as  desired,  an  analysis 
based  on  that  of  Reference  [4]  is  given.  Reference  [4]  considers  three  alternative  longitudinal  and 
lateral/directional  limited  authority  back-up  control  systems  for  the  F-4C.  The  research  of 
Reference  [4]  was  undertaken  to  develop  analytical  techniques  for  determining  the  minimum 
levels  for  the  back-up  control  system  parameters.  A  companion  project  utilized  piloted,  fixed- 
base  simulation  to  empirically  determine  the  minimum  levels  for  the  back-up  control  system 
parameters.  A  PIO  analysis  is  conducted  in  the  reference  by  applying  sinusoidal  input  describing 
function  theory.  Susceptibility  to  initiating  an  oscillation  is  investigated  by  assuming  the  pilot 
controls  in  a  normal  way  except  the  amount  of  lead  equilization  is  reduced.  The  ability  to  sustain 
the  oscillation  is  investigated  assuming  the  pilot  behaves  in  a  synchronous  manner.  The  dynamic 
model  contains  two  types  of  nonlinear  elements:  actuator  rate  and  displacement  limits  for  a  fully 
powered  manual  control  system,  and  direct  manual  control  authority  limits  determined  by  pilot 
strength  or  by  control  surface  travel  stops. 
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Figure  4-13  shows  the  block  diagram  for  the  longitudinal  back-up  control  system  A.  The 
shorthand  notation  used  in  Figure  4-13  is  (ca)  =  (s+ca).  The  parameters  for  the  pilot  loop  closures 
for  the  configuration  in  Figure  4-13  are  given  in  Table  4-4.  The  pilot  model  is  used  as  given  in 
Reference  [4]  and  was  developed  using  techniques  based  on  the  crossover  pilot  model. 


Table  4-4.  Model  Parameters  for  the  F-4C. 


Ky  =30.0  unit  thrust/(ft/sec) 
Kd=-0.03  deg/ft 
Xg  =0.2  sec 
(O3  =  1 0  rad/sec 


Ke  =-0.393  Kj  =-0.244  sec 

Kj  =-0.0015  deg/(ft/sec)  Xj  =0.33  sec 
(o„  =10  rad/sec  toj  =0.5  rad/sec 

(O4  =2.25  rad/sec 


Figure  4-14  shows  the  deflection  and  rate  limited  actuator  model  used  in  the  describing 
function  analysis.  Reference  [4]  develops  a  single  SIDF  for  the  combination  of  the  two  nonlinear 
elements  represented  in  Figure  4-14  and  uses  graphical  means  to  predict  the  existence  of  limit 
cycles.  However,  the  methods  of  this  section  can  be  employed  to  treat  the  two  nonlinear  elements 
as  separate  entities.  The  SIMULINK™  diagram  needed  to  define  M(s)  is  given  in  Figure  4-15. 

The  first  inport/outport  pair  corresponds  to  the  actuator  rate  limit  while  the  second  inport/outport 
pair  corresponds  to  the  actuator  deflection  limit. 

The  results  of  applying  the  limit  cycle  search  algorithms  given  in  Section  4. 3  .4. 3  are  presented 
in  Tables  4-5  and  4-6.  Also  shown  in  the  tables  are  the  results  from  the  given  reference.  A  dash 
in  a  column  indicates  that  no  limit  cycle  is  predicted  for  the  corresponding  configuration.  The 
tables  show  results  for  all  configurations  studied  in  Reference  [4]  though  only  longitudinal  back¬ 
up  control  system  A  has  been  presented  in  any  detail  here.  The  results  found  with  the  current 
method  agree  extremely  well  with  those  given  in  the  reference.  The  newly  calculated  amplitudes 
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Figure  4-13.  Longitudinal  Block  Diagram  for  F-4C  Back-up  Control  System  A. 
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Figure  4-14.  Rate  and  Position  Limited  Actuator. 


vary  by  less  than  5%  from  those  given  except  for  lateral/direction  system  C  which  varied  by  1 0%. 
The  predicted  frequencies  vary  by  less  than  3%  from  those  given  in  the  reference.  Since  the 
amplitudes  have  been  scaled  such  that  values  greater  than  unity  indicate  limiting  of  the  saturation 
elements,  it  is  seen  that  only  lateral/directional  system  A  is  simultaneously  rate  and  position 
limited. 


Table  4-5.  Comparison  of  Limit  Cycle  Results  for  the  Longitudinal  Systems. 

Reference  [41 _  _ Multivariable  Analysis _ 

Back-up  Ts  Ke  Rate  Deflection  co  Rate  Deflection  o  Phase  Difference 

System  (see)  fsec)  Aitip. _ Amp. _ (rad/s)  Amp. _ Amp.  (rad/s) _ (deg) - 


A 

0.2 

0 

12.3 

- 

1.1 

12.2 

0.79 

1.07 

-90 

A 

0 

0 

10.0 

- 

1.2 

10.4 

0.71 

1.19 

-90 

B 

0.2 

0 

- 

1.3 

1.4 

- 

1.30 

1.43 

- 

B 

0 

0 

- 

- 

- 

- 

“ 

- 

Table  4-6. 

Comparison  of  Limit  Cycle  Results  for  the  Lateral/Directional  Systems. 

Reference  [41 

Multivariable  Analysis 

Back-up 

^a 

Ka 

Rate 

Deflection 

0) 

Rate 

Deflection 

ca 

Phase  Difference 

System 

(sec) 

Amp. 

Amo. 

frad/s) 

mm 

Amo. 

trad/s) 

_ (deg) _ 

A 

0.2 

0 

15.0 

1.0 

0.14 

14.3 

1.01 

0.14 

-4 

B 

0.2 

0 

- 

4.3 

0.13 

- 

4.4 

0.13 

- 

C 

0.2 

0 

- 

5.0 

0.23 

- 

4.5 

0.23 

- 

Figure  4-16  shows  the  simulation  results  for  longitudinal  system  A  with  =  0,2.  The 
resulting  amplitudes,  frequency,  and  phase  are  seen  to  be  as  predicted  by  the  multivariable  analysis 
techniques.  Figures  4-17  and  4-1 8  show  the  altitude-altitude  rate  phase  plane  for  the  same 
configuration  as  Figure  4-16  for  small  and  large  initial  conditions  respectively.  The  resulting 
closed  curves  on  the  phase  plane  are  characteristic  of  limit  cycles  and  it  is  noted  that  the  same 
closed  curve  is  obtained  for  both  initial  conditions. 
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Input  Amplitudes  to  the  Nonlinear  Elements 


Figure  4-18.  Altitude-Altitude  Rate  Phase  Plane  (Large  Initial  Condition). 

While  Equation  (4-21)  can  be  searched  directly  using  a  nonlinear  equation  solver,  often  times 
a  good  initial  guess  is  necessary  to  ensure  convergence.  For  this  reason,  the  homotopy  method  as 
described  above  is  used  as  the  first  solution  technique.  Since  the  homotopy  algorithm  can  be  slow 
to  converge  when  using  numerical  partial  derivatives,  and  since  nonlinear  equation  solvers  can  be 
quite  rapid  when  in  the  neighborhood  of  a  solution,  a  combination  of  the  two  methods  can  be 
useful  in  obtaining  the  solution  in  the  following  manner.  An  initial  guess  is  obtained  from 
examining  the  stability  robustness  results.  The  homotopy  equations  are  integrated  using  a  low 
tolerance  to  yield  a  new  initial  guess.  The  nonlinear  equation  solver  can  then  be  employed  to 
yield  a  more  accurate  solution. 
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4.3.4.5  Summary 

Describing  function  analysis  provides  more  information  about  instability  than  structured 
singular  value  analysis  because  more  information  is  known  about  how  the  model  parameters  vary. 
The  structured  singular  value  calculation  yields  a  bound  on  the  largest  parameter  variation  or  the 
largest  sector  bound  of  a  nonlinearity.  The  describing  function  analysis,  on  the  other  hand,  yields 
the  magnitude  and  phase  of  resulting  oscillations.  However,  the  describing  function  analysis 
requires  a  fairly  sophisticated  search  for  solutions  whereas  computationally  efficient 
approximations  to  the  structured  singular  value  are  well-known.  As  a  result,  these  two  analysis 
methods  are  very  complementary  to  the  study  of  pilot-in-the-loop  oscillations.  Also,  the  existence 
of  hard  self  excitations  may  help  explain  why  certain  aircraft  exhibit  acceptable  flying  qualities  for 
some  pilots  and  yield  severe  PIOs  for  others.  Recall,  hard  self  excitation  describes  the  condition 
where  a  certain  amplitude  or  waveform  of  an  input  signal  is  necessary  to  excite  the  limit  cycle. 

4.3.5  Unified  Single-Loop  Analysis 

This  report  section  will  establish  ties  between  existing  single,  closed-loop  PIO  analysis 
methods  and  multivariable  stability  robustness  methods.  The  purpose  of  this  discussion  is  not  to 
refute  any  existing  methods  but  to  show  how  these  PIO  analysis  methods  are  closely  related  to 
concepts  in  stability  robustness  theory.  In  fact,  the  existing  methods  will  form  a  foundation  to 
establish  a  specific  new  single-loop  criteria. 

4.3. 5.1  Develop  Pilot  Model 

For  the  examples  of  this  section,  three  different  types  of  pilot  models  are  considered.  Both 
synchronous  and  MOCM  pilot  models  are  used  in  the  stability  robustness  analysis  of  Section 
4.3 .5.3 .  For  the  synchronous  (pure  gain)  pilot  model,  the  pilot  gains  were  chosen  such  that  the 
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pilot/vehicle  transfer  functions  have  at  least  6dB  of  gain  margin  and  45  degrees  of  phase  margin. 
The  values  for  the  MOCM  model  are  as  given  in  Table  4-1  and  were  chosen  as  in  Reference  [11] 
with  the  disturbance  filter  taken  as  l/(s+l)2.  Smith’s  nonlinear  contactor  pilot  model  is  used  for 
the  describing  function  analysis  of  Section  4. 3. 5. 4. 

4.3. 5.2  Identification  of  the  Dynamic  Elements  Suspected  of  Contributing  to  PIO 

PIO  susceptibility  will  be  examined  in  the  face  of  variations  in  the  pilot  model.  For  this 
section,  the  variation  will  be  taken  as  a  multiplicative  uncertainty  in  pilot  equalization. 

4.3. 5.3  Stability  Robustness  Analysis 

The  standard  form  M-A  diagram  is  shown  in  Figure  4-2.  A  feedback  loop  uncertainty  or 
variation  is  represented  by  the  matrix  A  while  M(s)  represents  the  linear,  lumped  parameter  part  of 
the  feedback  system.  The  harmonic  balance  equations  are  found  by  substituting  s=ja)  in  (4-10)  to 
yield, 

[I  +  M(jco)A]  xOco)  =  0  (4-25) 

Note  that  A  may  also  be  a  function  of  frequency,  i  .e.  A((o).  Since  we  are  interested  in  non-trivial 
solutions  (x(j(o)  ^  0),  one  can  search  for  solutions  of  (4-25)  above  from,['^k42] 

det[I  +  M0a))A]  =  0  (4-26) 

For  this  section  only,  we  will  consider  only  a  single  variation  parameter.  Thus,  the  matrix  A  is 
actually  a  scalar  variable  that  will  be  called  6ni  and  M(s)  is  a  single-input,  single-output  transfer 
function  that  will  be  represented  by  m(s).  Under  these  conditions.  Equation  (4-26)  becomes, 

1  +  m(ja))6m  =  0  (4-27) 

or, 

m(j(D)  =  -l/Sm  (4-28) 
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Now  consider  the  block  diagram  shown  in  Figure  4-19.  This  block  diagram  represents  a 
pilot/aircraft  closed-loop  system.  The  transfer  function  g(s)  models  the  aircraft  while  po(s) 
models  the  (nominal)  pilot.  The  parameter  6  will  represent  in-line  multiplicative  variations  in  the 
pilot  equalization.  These  variations  might  come  from  inappropriate  pilot  action,  for  example. 
One  might  think  of  the  actual  pilot  p(s)  as  a  multiplicative  combination  of  a  nominal  model  po(s) 
and  a  variation  6,  i.e.  p(s)  =  6po(s). 


pilot  variation  aircraft 


Figure  4-19.  Single  Loop  Pilot/Aircraft  System. 

The  governing  equations  for  the  block  diagram  of  Figure  4-19  are, 

y(s)  =  5x(s) 

x(s)  =  -po(s)g(s)y(s) 

Substituting  (4-29)  into  (4-30)  yields, 

[1  +  po(s)g(s)5]x(s)=  0 

so  that  the  harmonic  balance  equations  for  this  case  are, 

[  1  +  po0®)g(ja))6]x(j(o)  =  0 

Ignoring  the  trivial  solution  of  (4-32),  we  are  therefore  interested  in  solutions  to, 

1  +  po(jo))g(ja))5  =  0 

or, 


(4-29) 

(4-30) 

(4-31) 

(4-32) 

(4-33) 
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Po(ja>)gO®)  =  -1/5  (4-34) 

Comparing  (4-34)  and  (4-28),  one  should  quickly  see  that  m(s)  is  simply  the  open-loop 
transfer  function  Po(s)g(s).  The  pilot  model  variation  parameter  6  is  identical  to  5ni.  Equation  (4- 
34)  can  be  used  to  describe  several  existing  PIO  analysis  methods. 

Let  the  variation  parameter  in  Figure  4-19  be  an  arbitrary  change  in  pilot  gain.  Let  5  be 
written  as  6  =  K5,  where  K5  represents  the  change  in  pilot  gain.  Assuming  a  synchronous  nominal 
pilot  model  Kp,  the  harmonic  balance  equation  (4-34)  is  then, 

KpgOco)  =  -1/K8  (4-35) 

Unfortunately,  we  do  not  have  a  specific  value  of  K5  to  test  (4-35)  for  a  given  airplane.  It  is 
expected  that  the  actual  pilot  gain  variation  would  occur  when  the  pilot  becomes  excited  or 
agitated  during  a  PIO  event.  One  might  expect  that  an  aircraft  that  resists  changes  in  pilot  gain 
might  be  less  susceptible  to  PIO  than  an  airplane  that  is  very  sensitive  to  pilot  gain  changes. 

The  stability  test  implied  by  (4-35)  is  simply  a  test  for  adequate  gain  margin  of  the 
pilot/aircraft  system.  This  test  can  be  performed  using  a  Nyquist  plot  or  a  Bode  diagram  as 
shown  in  Figure  4-20  for  Configurations  H2-1  and  H2-5  of  the  HAVE  PIO  database.  Using  the 
synchronous  pilot  model,  the  nominal  pilot  gains  for  Configurations  H2-1  and  H2-5  were  chosen 
as  Kp  =  1.24  and  Kp  =  1.09  respectively. 

The  first  thing  to  notice  in  Figure  4-20  is  that  the  gain  margin  of  Configuration  H2-1  is  much 
larger  than  H2-5.  The  average  PIO  rating  for  Configuration  H2-1  was  1 .0  while  the  average 
rating  for  H2-5  was  4.33  (a  rating  above  2.0  is  considered  PIO  susceptible).  Consequently,  one 
might  claim  that  the  gain  margin  is  a  good  indicator  of  PIO  susceptibility.  Figure  4-21  shows  the 
gain  margin  for  each  of  the  eighteen  HAVE  PIO  configurations  plotted  with  the  average  PIO 
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Figure  4-20.  Bode  Plots  for  HAVE  PIO  Configurations  H2-1  and  H2-5. 

rating  given  during  flight.  As  expected.  Figure  4-21  shows  that  the  PIO  rating  has  a  definite 
correlation  with  gain  margin.  The  configurations  with  large  relative  gain  margins  tend  to  have 
lower  PIO  ratings.  However,  one  can  also  see  from  Figure  4-21  that  some  configurations  with 
small  relative  gain  margins  (on  the  order  of  2-3  absolute)  still  yield  low  PIO  ratings. 

Figure  4-20  also  shows  that  the  phase  margins  for  Configuration  H2-1  and  H2-5  are  nearly 
identical.  The  identical  phase  margins  stem  from  the  criteria  used  to  choose  the  nominal 
synchronous  model  pilot  gain  Kp  (to  insure  at  least  45  deg  phase  margin  and  6  dB  gain  margin). 
Thus,  using  this  method  for  choosing  pilot  gain  leads  to  the  conclusion  that  phase  margin  does  not 
seem  to  correlate  with  PIO  rating. 


50 


7 


Figure  4-21.  HAVE  PIO  Configuration  Gain  Margins. 


Hess  and  Kalteis  have  developed  a  PIO  criterion  that  explicitly  defines  the  required  magnitude 
slope  of  the  pilot/aircraft  open-loop  system,  The  recommended  pilot  model  is  an  Optimal 
Control  Model  (OCM).  The  Hess-Kalteis  Criterion  is  shown  graphically  in  Figure  4-22.  Note 
that  the  boundary  limits  the  minimum  gain  crossover  frequency  to  2.0  rad/s.  For  a  pilot/aircraft 
open-loop  response  with  a  crossover  of  exactly  2.0  rad/s,  the  criterion  also  limits  the  magnitude 
slope  to  at  least  -20  dB/decade  in  the  crossover  frequency  range.  Thus,  like  the  open-loop  criteria 
of  Smith-Geddes  and  Hoh,  the  Hess-Kalteis  criterion  tends  to  identify  configurations  with  steeper 
magnitude  slopes  in  the  crossover  frequency  range  as  PIO  resistant. 

The  military  flying  qualities  specification  includes  the  Neal-Smith  closed-loop  criteria  for 
analysis  of  longitudinal  short-period  dynamics,  l^^l  A  complete  lead/lag  (crossover)  pilot  model  is 
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FREQUENCY  RAO/SEC 


Figure  4-22.  Hess-Kalteis  PIO  Criterion. 

specified,  and  measures  such  as  bandwidth  and  closed-loop  resonant  peak  are  used  to  determine 
the  aircraft's  flying  qualities  rating. 

A  variation  of  the  Neal-Smith  criteria  can  also  be  developed  using  the  M-A  standard  form. 
Using  the  harmonic  equations  from  (4-33),  add  and  subtract  PoO®)gO®)  to  yield, 

1  +  po(jo))g0®)5  +  PoO®)g(j®)  -  PoOo>)g(jfl>)  =  0  (4-36) 

Now  combine  terms  in  (4-36)  so  that, 

[1  +  po(j®)gGtt')]  +  (5  -  I)po0to)g0w)  =  0  (4-37) 

and  dividing  by  [1  +  poO<^)gO®®)]  reveals. 


(5-l)po0cQ)g(i(D) 

[1+Po(jco)g0a))] 


(4-38) 


or, 

Po(jco)g(io)) 

[l+Po(jcf))g(jca)]  5-1 

From  (4-39),  one  should  see  that  now  Sm  =  1/6-1 .  The  transfer  function  m(s)  is  equal  to 
Po(s)g(s)/(l+po(s)g(s)),  which  is  the  same  as  the  closed-loop  transfer  function  from  rc(s)  to  r(s)  in 
Figure  4-19. 
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The  Neal-Smith  flying  qualities  criteria  specifies  the  maximum  allowable  resonant  peak  of  the 
closed-loop  frequency  response.  For  Level  1  airplanes,  the  resonant  peak  must  be  less  than  3  dB 
and  for  Level  2,  the  resonant  peak  must  be  less  than  9  dB.  Let  this  resonant  peak  be  denoted  by 


P,  therefore, 

max  Po(io>)gO«>) 

0)  ll+po(j®)gOo>)l 

If  5  is  not  a  function  of  frequency,  (4-39)  is  satisfied  whenever. 


(4-40) 


(4-41) 


Now  suppose  5  represents  only  a  change  in  pilot  gain,  i.e.  5  =  Ks  as  before.  When  5  =  Kf>,  (4-41 ) 
can  be  solved  for  K5  in  terms  of  P  to  yield, 

K6  =  1  +  1/P  (4-42) 

Using  (4-42),  the  Neal-Smith  criteria  can  be  interpreted  as  a  limitation  on  the  maximum 
tolerable  increase  in  pilot  gain.  For  Level  1  aircraft,  (4-42)  reveals  that  when  P  =  3  dB,  then  Ks  = 
1.71  (4  65  dB).  This  result  means  that  the  Level  1  Neal-Smith  criteria  allows  up  to  a  1.71  factor 
increase  in  pilot  gain.  Not  surprisingly,  the  Level  2  criteria  of  P  =  9  dB  translates  to  a  maximum 
pilot  gain  increase  of  Ks  =  1.35  (2.64  dB).  Thus,  a  Level  2  airplane  can  therefore  tolerate  less 

pilot  gain  change  than  a  Level  I  airplane. 

In  summary,  each  of  the  existing  single,  open-loop  and  closed-loop  criterion  appear  to  focus 
upon  various  measures  of  frequency  response  shape.  It  is  also  important  to  note  that  each  of 
these  existing  criteria  can  be  motivated  by  how  arbitrary  pilot  gain  changes  can  affect  closed-loop 
stability.  The  Smith-Geddes  and  Hess-Kalteis  methods  deal  mainly  with  the  pilot/aircraft 
magnitude  slope  near  the  actual  (Hess-Kalteis)  or  implied  (Smith-Geddes)  gain  crossover 
frequency.  Changes  in  pilot  gain  will  affect  the  pilot/aircraft  gain  crossover  frequency  directly. 
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The  Neal-Smith  closed-loop  criterion  can  also  be  interpreted  as  limiting  the  maximum  pilot  gain 
variation.  Therefore,  each  of  these  methods  would  appear  to  relate  directly  to  stability  robustness 
with  respect  to  changes  or  variations  in  pilot  gain. 

Pilot  gain  also  affects  phase  margin.  Since  phase  margin  is  defined  relative  to  the  -180  degree 
point,  it  is  reasonable  to  develop  criteria  governing  the  shape  of  the  aircraft  phase  response  near 
the  phase  angle  of -180  deg.  Hoh's  bandwidth  criteria  directly  specifies  a  minimum  phase  slope 
near  (0180  while  the  Smith-Geddes  criteria  indirectly  favors  shallow  phase  slopes. 

It  has  already  been  shown  that  gain  margin  correlates  fairly  well  with  PIO  rating  for  the 
HAVE  PIO  configurations.  Since  most  of  the  existing  criteria  studied  herein  have  direct  ties  to 
gain  margin  and,  to  a  lesser  extent,  phase  margin  concepts,  it  seems  worthwhile  to  consider 
methods  which  can  measure  simultaneous  gain  and  phase  variations.  The  idea  of  a  vector  stability 
margin  is  therefore  offered  as  a  unified  measure  of  PIO  susceptibility.  The  vector  margin  was 
originally  developed  as  a  stability  robustness  measure  for  linear  systems, 

Figure  4-23  shows  the  graphical  definition  of  a  vector  stability  margin  for  the  pilot/vehicle 
representation  shown  in  Figure  4-19.  The  vector  margin  P  is  defined  as  the  minimum  distance 
from  the  graph  of  poO(o)g(jw)  to  the  critical  -1+jO  point.  One  can  see  that  p  is  related  to 
traditional  gain  and  phase  margins.  Note  from  Figure  4-23  that  the  vector  p  is  smaller  then  Pi, 
which  is  the  vector  which  defines  the  traditional  gain  margin.  It  is  also  smaller  than  P2,  which  is 
the  vector  which  defines  the  traditional  phase  margin.  Therefore,  when  P  =  pi,  the  length  of  the 
vector  P  can  be  used  to  give  an  exact  representation  for  traditional  gain  margin.  Conversely, 
when  P  =  P2,  the  length  of  P  can  be  used  to  find  the  exact  value  of  the  traditional  phase  margin 

(e  g.  PM  =  2sin-HP2/2)). 
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Figure  4-23.  Vector  Margin  Definition. 

The  vector  stability  margin  has  many  appealing  characteristics  when  considered  for  PIO 
analysis.  The  vector  margin  deals  with  the  shape  of  the  pilot/aircraft  frequency  response  near  the 
critical  point  in  a  manner  very  similar  to  the  various  existing  PIO  criteria.  In  addition,  the  vector 
margin  automatically  establishes  a  critical  frequency  as  the  frequency  to  where  Po0^o)g(jt>>)  is 
nearest  the  -1+jO  point.  Note  that  this  frequency  point  is  not  necessarily  the  open-loop  gain 
crossover  frequency  or  the  frequency  point  associated  with  -180  phase  angle.  It  is  hypothesized 
that  this  frequency  point  will  be  closely  related  to  the  frequency  of  PIO. 

The  vector  margin  can  be  readily  developed  from  the  M-A  concepts  described  previously. 
Consider  (4-33)  again,  add  and  subtract  5  to  yield, 

1  +  po(jco)g(ja))6  +  6-6  =  0  (4-43) 

Combine  terms  in  (4-43), 

( 1  -  6)  +  [  1  +  po0w)g0o>)]S  =  0  (4-44) 
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and  divide  by  6, 


[1  +  poOffl)gO®)l  =  ^  =  1-1/5  (4-45) 

Comparing  (4-45)  to  (4-28),  we  note  that  m(s)  =  1  +  po(s)g(s)  and  6ni  =  6-1/5. 

The  vector  margin  is  a  vector  pointing  from  po(jo))g(j(o)  to  the  -1+jO  point.  Since  we  are 
interested  in  the  smallest  distance,  let  P  be  the  smallest  magnitude  of  [1  +  po0®)80®)]  or, 

P  =  |1  +  Po(jco)gOco)|  (4-46) 

We  will  also  define  the  critical  frequency  as  the  frequency  where  the  minimum  is  obtained 
in  (4-46). 

The  value  of  P  is  shown  graphically  in  Figure  4-23.  The  aircraft  depicted  in  Figure  4-23  is  the 
B-2  pitch  attitude  response  during  air-to-air  refueling,  Using  a  synchronous  pilot  model  gain 
of  Kp  =  20,  the  length  of  the  vector  margin  in  this  case  is  p  =  0.457  at  a  frequency  of  8  . 13  rad/s. 
For  reference,  this  aircraft  was  known  to  experience  PIOs  in  the  configuration  under  study. 

Figure  4-24  shows  a  magnitude  plot  of  1  +  po(jci))g(j(o)  for  HAVE  PIO  Configurations  H2-1 
and  H2-5.  The  nominal  pilot  model  is  again  chosen  to  be  a  simple  gain  po(s)  =  Kp.  The  value  of 
P  will  be  the  minimum  point  on  this  frequency  response.  From  Figure  4-24,  we  can  see  that  P  = 
0.54  for  Configuration  H2-1  and  P  =  0.43  for  Configuration  H2-5.  The  critical  frequencies  are 
4. 1  and  2.0  rad/s  for  Configurations  H2-1  and  H2-5,  respectively. 

As  expected,  the  value  of  p  is  smaller  for  H2-5  than  H2-1 .  This  characteristic  means  the 
Nyquist  plot  of  poG®“)gO<^)  comes  closer  to  the  -1+jO  point  for  Configuration  H2-5  than  for  H2-1 . 
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Frequency  (rad/sec) 


Figure  4-24.  Vector  Margins  for  H2-1  and  H2-5  (Synchronous  Pilot). 

It  also  means  that  Configuration  H2-1  can  sustain  larger  simultaneous  pilot  gain  and  phase 
changes  than  can  H2-5. 

Table  4-7  summarizes  the  results  of  computing  the  vector  margin  for  each  of  the  HAVE  PIO 
Configurations.  A  nominal  synchronous  gain  pilot  model  was  used  for  these  results.  Again,  the 
pilot  model  gain  Kp  was  chosen  as  the  minimum  value  which  leads  to  at  least  6  dB  gain  margin 
and  at  least  45  deg  phase  margin.  Also  shown  in  Table  4-7  is  the  critical  frequency  where  P  is 
found. 

Figure  4-25  shows  the  vector  margin  plotted  with  the  average  PIO  rating  given  in  flight  for 
each  HAVE  PIO  configuration.  The  trend  of  increasing  PIOR  with  decreasing  vector  margin  is 
clearly  evident  in  Figure  4-25. 
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Vector  Margin  (^) 


Table  4-7.  Vector  Margin  using  Synchronous  Pilot. 


Configuration 

PIOR 

Kp 

P 

“crit  (rad/s) 

H2-B 

2.0 

1.88 

0.57 

7.15 

H2-1 

1.0 

1.24 

0.54 

4.10 

H2-5 

4.33 

1.09 

0.43 

2.05 

H2-7 

3.0 

0.93 

0.43 

3.33 

H2-8 

4.0 

0.78 

0.45 

3.11 

H3-D 

1.0 

1.83 

0.49 

6.67 

H3-1 

2.33 

2.20 

0.51 

8.50 

H3-3 

1.67 

1.38 

0.45 

4.25 

H3-6 

2.0 

1.27 

0.47 

5.81 

H3-8 

3.67 

0.95 

0.48 

4.71 

H3-12 

4.5 

0.67 

1.98 

H3-I3 

4.5 

0.70 

0.46 

2.61 

H4-1 

1.0 

1.54 

5.41 

H4-2 

1.33 

1.16 

0.45 

4.25 

H5-1 

1.0 

1.27 

0.55 

2.90 

H5-9 

4.0 

0.81 

2.12 

H5-10 

5.0 

0.62 

0.45 

1.91 

H5-I1 

3.0 

0.91 

0.45 

2.35 

0  1  2  3  4  5  6 


Flight  Test  PIOR 


Figure  4-25.  Vector  Margin  and  PIO  Rating  (Synchronous  Pilot). 
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Figure  4-26  shows  the  critical  frequency  and  the  actual  PIO  frequency  for  the  thirteen  HAVE 
PIO  configurations  that  actually  resulted  in  a  PIO  event.  The  correlation  of  the  vector  margin 
critical  frequency  and  the  actual  PIO  frequency  is  good  These  results  show  that  the  actual  PIO 
frequency  is  slightly  higher  than  the  critical  frequency  predicted  by  the  vector  margin. 

Figure  4-27  shows  the  critical  frequency  plotted  against  the  vector  margin  using  the 
synchronous  pilot  model.  Note  that  with  only  two  exceptions  (H3-1  and  H3-8)  the  PlO-prone 
configurations  (PIOR  >  2.0)  are  grouped  together  with  small  vector  margins  and  low  critical 
frequencies.  The  PlO-free  configurations  tend  to  have  larger  vector  margins  and  critical 
frequencies. 


2  4  6  8  10  12 

Critical  Frequency  SP  (rad/s) 

Figure  4-26.  Predicted  PIO  Frequencies  (Synchronous  Pilot). 
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Critical  Frequency  SP  (rad/s) 

Figure  4-27.  Vector  Margin  and  Critical  Frequency  (Synchronous  Pilot). 

The  two  exceptions  to  the  trends  in  Figure  4-27  are  configurations  H3-1  and  H3-8. 
Configuration  H3-1  has  a  short-period  natural  frequency  of  4.1  rad/s  and  damping  ratio  of  1.0. 
No  other  dynamics  are  included  in  H3-1,  yet  the  configuration  was  flown  three  times  with  ratings 
of  3,  2,  and  2.  The  PIO  frequency  experienced  during  flight  test  was  10.4  rad/s.  As  a  result,  it 
would  appear  that  this  configuration  may  have  encountered  a  higher-fi-equency  phenomenon  than 
is  considered  herein.  Other  PIO  criteria  have  also  had  difficulty  predicting  the  behavior  of 
H3-I.[n] 

A  nominal  synchronous  pilot  model  was  used  in  the  preceding  results.  The  vector  margin  can 
also  be  computed  using  more  sophisticated  pilot  models.  Table  4-8  shows  the  results  of  using  the 
more  sophisticated  MOCM  operator  model  described  in  Section  4.3.1 
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Table  4-8.  Vector  Margin  using  MOCM  Pilot. 


Configuration  PIOR 


H2-B 

2.0 

H2-1 

1.0 

H2-5 

4.33 

H2-7 

3.0 

H2-8 

4.0 

H3-D 

1.0 

H3-1 

2.33 

H3-3 

1,67 

H3-6 

2.0 

H3-8 

3.67 

H3-12 

4.5 

H3-13 

4.5 

H4-1 

1.0 

H4-2 

1.33 

H5-1 

1.0 

H5-9 

4.0 

H5-10 

5.0 

H5-11 

3.0 

P 

cocrit  (rad/s) 

0.433 

6.44 

0.429 

5.81 

0,408 

4.40 

0.412 

5.23 

0.410 

5.05 

0.423 

6.44 

0,429 

6.67 

0.414 

5.61 

0.415 

6.01 

0.410 

5.61 

0.405 

4.25 

0.408 

4.55 

0.425 

6.22 

0.413 

5.61 

0.436 

5.42 

0.408 

4.55 

0.405 

4.25 

0.413 

4.71 

Notice  that  the  values  for  p  are  much  smaller  using  the  MOCM  operator  model  as  opposed  to 
the  synchronous  model.  This  result  occurs  because  the  MOCM  model  includes  neuromuscular 
dynamics  as  well  as  time  delay.  These  pilot  limitations  have  the  effect  of  reducing  the  vector 
stability  margins.  In  a  sense,  the  synchronous  pilot  model  leads  to  a  somewhat  optimistic 
prediction  of  pilot/aircraft  stability  robustness. 

Figure  4-28  shows  the  MOCM  vector  margins  plotted  with  their  associated  critical 
frequencies  for  the  HAVE  PIO  configurations.  Figure  4-28  now  shows  that  the  PlO-prone 
configurations  are  all  grouped  together  with  only  the  exception  of  H3-1 .  These  configurations  are 
grouped  in  an  area  where  the  vector  margin  is  less  than  p  s  0.415  and  the  critical  frequency  is  less 
than  approximately  5.5  rad/s. 
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Figure  4-28.  Vector  Margin  and  Critical  Frequency  (MOCM  Pilot). 

To  help  establish  the  value  of  P  that  defines  PlO-resistant  and  PlO-prone  aircraft,  one  can 
consider  the  relationship  of  the  vector  margin  to  the  Neal-Smith  criteria.  If  only  pilot  gain 
variations  are  considered  so  5  =  Kg,  then  (4-45)  and  (4-46)  become, 

P  =  1  -  I/K5  (4-47) 

It  was  previously  found  that  the  Neal-Smith  criteria  would  allow  a  pilot  gain  variation  of  up  to 
Kg  =  1.71  for  Level  1  flying  qualities.  Using  equation  (4-47),  we  find  that  the  Neal-Smith  Level  1 
criteria  yields  the  value  of  P  =  0.41 5. 

Figure  4-29  shows  the  critical  frequency  predicted  by  the  vector  margin  and  the  MOCM 
plotted  along  with  the  actual  PIO  frequencies.  The  correlation  in  frequency  is  fairly  good; 
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Figure  4-29.  Predicted  PIO  Frequencies  (MOCM  Pilot  Model). 


however,  the  critical  frequencies  are  predicted  too  large  when  the  actual  value  is  small  and  vice 
versa.  Nevertheless,  the  trend  of  predicted  critical  frequency  and  actual  PIO  frequency  is  clearly 
evident. 

4.3. 5.4  Multivariable  Describing  Function  Analysis 

The  preceding  analysis  has  studied  the  effect  of  pilot  gain  and  phase  variations  on  closed-loop 

stability.  One  particular  pilot  gain  variation  occurs  when  the  pilot  response  is  limited.  Now, 
consider  the  case  where  the  pilot  is  to  be  modeled  by  Smith’s  contactor  pilot  model. l^^l  Since  we 
are  searching  for  instances  of  sinusoidal  behavior,  the  contactor  can  be  modeled  using  its 
sinusoidal  describing  function  in  Table  4-3, 

N(a,(D)=”  (4-48) 
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where  it  has  been  assumed  that  the  input  to  the  contactor  is  x(t)  =  asin((ot)  and  the  contactor 
output  level  is  'a'.  Note  that  whenever  the  input  amplitude  is  smaller  than  the  output  level,  i.e. 
a/a  >  1 ,  the  describing  function  magnitude  is  greater  than  unity.  In  this  case,  the  contactor  output 
level  might  represent  the  maximum  control  stick  deflection  or  force  level  available  to  the  pilot. 

Returning  to  Figure  4-19,  we  will  assume  that  the  nominal  pilot  is  modeled  by  a  simple  gain 
Kp  (a  synchronous  pilot  model).  The  contactor  will  represent  the  variable  element  5.  As  a  result, 
the  stability  condition  (4-34)  becomes, 

PoOco)gO®)  = -1/5  =>  Kpg(jcD)  =  -l/N(a,(o)  =  -”  (4-49) 

The  existence  of  solutions  to  (4-49)  can  be  found  using  a  Nyquist  diagram.  Figure  4-30  shows  a 
Nyquist  diagram  of  Kpg(jo))  where  g(jo))  models  Configuration  H2-5  of  the  HAVE  PIO 

database.t^l  The  nominal  pilot  gain  is  as  previously  chosen,  Kp  =  1.09. 

The  Nyquist  plot  in  Figure  4-30  shows  the  two  sides  of  (4-49),  i.e.  Kpg(jto)  and  -l/N(a,(D). 
The  plot  of -1/N(a,c())  assumes  that  the  ratio  a/a  varies  from  0.2  to  2.0.  When  the  two  lines  in 
Figure  4-30  cross,  a  solution  to  (4-49)  is  obtained.  Figure  4-30  shows  that  (4-49)  has  a  solution 
when  CO  =  2.4  rad/s  and  ot/a  =  0.56.  This  solution  corresponds  to  a  stable  limit  cycle  of  the 
pilot/aircraft  closed-loop  system.  Figure  4-3 1  shows  the  phase  plane  trajectory  for  this  limit  cycle 
when  the  contactor  output  level  is  a  =  1 .  Note  that  the  response  starts  with  a  very  small  initial 
condition  and  builds  to  a  constant  magnitude  oscillation  which  can  be  likened  to  a  fiilly  developed 
PIO.  The  contactor  describing  function  in  (4-49)  amounts  to  a  change  in  pilot  gain  with  respect 
to  its  nominal  value.  The  gain  change  has  a  specific  relationship  to  input  magnitude  and  can  lead 
to  possible  oscillatory  behavior. 
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Figure  4-30.  Nyquist  Plot  for  Contactor  Pilot  Model 


4. 3. 5. 5  Summary 

Existing  PIO  criterion  seem  to  focus  on  various  measures  of  the  frequency  response  shape. 
Furthermore,  each  of  these  existing  criteria  can  be  motivated  by  how  arbitrary  pilot  dynamic 
changes  can  affect  closed-loop  stability.  Therefore,  each  of  these  methods  would  appear  to  relate 
directly  to  stability  robustness  with  respect  to  changes  or  variations  in  pilot  gain.  The  vector 
stability  margin  is  offered  as  one  unifying  measure  of  the  shape  of  the  pilot/aircraft  frequency 
response  near  the  critical  point  and  as  an  indication  of  PIO  susceptibility. 

4.3.6  Unified  Multi-Loop  Analysis 

The  vector  margin  concept  is  readily  extended  to  multi-loop  systems  because  the  single-loop 
margin  previously  discussed  is  actually  a  special  case  of  a  more  general  theory. 

4.3.6. 1  Develop  Pilot  Model 

The  methods  of  pilot  modeling  given  in  the  previous  section  can  be  extended  for  use  in  the 
multi-loop  case.  However,  for  this  section,  only  a  generic  transfer  function  representation  of  the 
pilot  is  considered  without  reference  to  how  the  model  is  to  be  derived. 

4.3. 6. 2  Identification  of  the  Dynamic  Elements  Suspected  of  Contributing  to  PIO 

The  impact  of  simultaneous  gain  and  phase  variations  in  each  loop  of  the  pilot  model  will  be 

considered  as  the  primary  uncertainties  in  this  section.  No  specific  nonlinearities  are  considered. 

4.3. 6.3  Stability  Robustness  Analysis 

Consider  the  block  diagram  in  Figure  4-32.  The  transfer  function  matrix  G(s)  in  Figure  4-32 
represents  the  aircraft  transfer  functions.  Po(s)  is  a  matrix  of  transfer  functions  representing  the 
pilot.  The  matrix  A  represents  changes  in  pilot  equalization  and  is  assumed  to  be  diagonal.  This 
block  diagram  assumes  single,  individual  pilot  variations  in  each  pilot  manipulator  channel.  For 
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Figure  4-32,  Multi-Loop  Pilot/Aircraft  System. 

example,  a  lateral-directional  model  might  include  pilot  gain  and  phase  variations  in  the  lateral 
control  stick  and  rudder  pedal  channels. 

From  Figure  4-32,  the  loop  equations  are  (rc(s)  =  0), 

y(s)  =  Ax(s)  (4-50) 

x(s)  =  -Po(s)G(s)y(s)  (4-5 1 ) 

Combining  the  loop  equations  and  substituting  s=j{o  leads  to  the  expression, 

[I  +  Po(jco)GOa))A]  x(j(D)  =  0  (4-52) 

Comparing  (4-52)  to  (4-25)  shows  that  M(s)  is  equal  to  the  open-loop  return  ratio  matrix 
Po(s)G(s)  for  this  case.  Since  we  are  interested  in  non-trivial  solutions  to  (4-52),  we  seek  to  find 
solutions  of, 

I  +  Po(ja))GOco)A  =  0  (4-53) 

As  before,  add  and  subtract  A, 

(I  -  A)  +  [I  +  PoGto)G(jco)]A  =  0  (4-54) 

and  multiply  by  A’^  to  yield, 

[I  +  Po(j(o)GO(o)]  =  (A-I)A-l  (4-55) 
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Equation  (4-55)  is  essentially  the  multi-loop  equivalent  to  (4-45).  However,  since  (4-55)  is  a 
matrix  expression,  some  measure  of  matrix  norm  is  required  to  define  the  distance  to  the  critical 
-1+jO  point.  Most  often  the  matrix  singular  value  is  used  as  the  measure,  in  which  case  the  multi¬ 
loop  vector  margin  is  defined  as, 

P  =  Omin[I  +  PoO®)GO®)]  (4-56) 

where  a„,i„  is  the  minimum  matrix  singular  value. 

Equation  (4-56)  would  be  utilized  in  the  same  manner  as  the  single-loop  equivalent  in  (4-46). 
A  value  of  P  is  computed  from  (4-56)  and  checked  against  a  minimum  acceptable  value  (yet  to  be 
determined).  The  frequency  at  which  the  minimum  in  (4-56)  is  found  will  determine  the  predicted 
PIO  frequency. 

The  vector  margin  defined  by  (4-56)  has  been  used  for  several  years  to  define  the  gain  and 
phase  stability  margins  for  multi-loop  feedback  systems.  147]  Pq^  example,  let  A  =  diag(Kiei*f'i) 
where  i  is  the  loop  index.  The  right  hand  side  of  (4-55)  can  then  be  written  as, 

(A-I)A-l  =  I-A-1  =  diag(l-l/Kie)‘l>i)  =  diag(l-e-J‘l>i/Ki)  (4-57) 

Now  substitute  e*J‘t’  =  cos(|)  -  jsin(|)  into  (4-57)  and  find  the  largest  magnitude, 

|(l-cos(|)i/Ki)+jsin(l)i/Ki|  (4-58) 

After  further  manipulation,  one  finds  that  (4-55)  will  not  have  a  solution  as  long  as, 

P  >  Jil-UK,y  +  2(l-cos(t)i)/K,  (4-59) 

i 

Equation  (4-59)  establishes  a  relation  between  pilot  gain  variations  K5,  phase  variations  (j),  and 
the  vector  margin  p.  Any  combinations  of  K5  and  (])  that  satisfy  (4-59)  will  insure  that  the 
harmonic  conditions  in  (4-55)  are  not  met.  Therefore,  (4-59)  has  been  used  to  define  the  largest 
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gain  and  phase  margins  that  a  multi-loop  feedback  system  can  tolerate  before  a  stability  condition 
is  violated.  12<>]  For  application  to  PIO  analysis,  the  multi-loop  vector  margin  measures  the  largest 
simultaneous  gain  and  phase  variations  that  the  pilot  can  introduce  in  each  loop  channel  before 
oscillatory  behavior  of  the  closed-loop  system  is  possible. 

A  minimum  acceptable  value  of  P  for  the  multi-loop  case  is  not  yet  known.  It  is  likely  that 
PIO  resistant  aircraft  will  have  smaller  values  of  p  in  the  multi-loop  case  than  in  the  single-loop 
case  simply  because  several  pilot  loops  are  involved  at  once.  However,  the  minimum  singular 
value  measure  is  also  known  to  be  conservative  with  respect  to  unit  scaling.  Therefore,  a  scaling 
matrix  like  that  used  in  structured  singular  value  analysis  will  be  required.  These  subjects  will  be 
considered  in  future  research. 

4.3. 6. 4  Multivariable  Describing  Function  Analysis 

Once  the  nonlinear  elements  of  the  multi-loop  system  are  isolated,  a  describing  function 

analysis  can  be  carried  out  exactly  as  described  in  Section  4.3.4. 

4.3.6. 5  Summary 

The  vector  stability  margin  is  given  for  the  multi-loop  case.  It  is  suggested  that  the  vector 
margin  could  be  used  as  an  indication  of  PIO  susceptibility  in  a  manner  similar  to  that  of  Section 
4.3.5.  However,  a  minimum  acceptable  value  for  this  margin  is  not  yet  known. 


69 


4.3.7  Unifled  Multivariable  Analysis 

This  section  describes  a  unified  method  for  multivariable  PIO  analysis.  The  term  multivariable 
is  used  here  to  imply  that  the  variable  dynamics  might  appear  anywhere  in  the  pilot/vehicle 
feedback  system.  In  order  to  exercise  the  stated  multivariable  analysis  method,  it  was  desired  to 
use  the  techniques  to  analyze  an  aircraft  that  had  experienced  in  flight  PIOs.  For  the  analysis  to 
be  meaningful,  sufficient  modeling  data  needed  to  be  available  for  the  vehicle  along  with  flight 
time  histories  of  the  oscillation  events.  The  M2-F2  lifting  body  vehicle  was  chosen  as  an  example 
case  study.  All  information  for  the  M2-F2  was  taken  from  Reference  [48]. 

The  M2-F2  was  a  half  cone  lifting  body  designed  as  part  of  the  manned  lifting  body  flight  test 
program  conducted  in  the  1 960s.  Basic  research  was  being  conducted  on  manned  reentry  vehicles 
designed  to  maneuver  and  land  unpowered  at  a  specified  location.  The  M2-F2  had  a  conventional 
fighter  aircraft  type  cockpit  with  aerodynamic  control  being  provided  by  upper  eleven  flaps,  a 
lower  flap,  and  twin  rudders.  A  three  view  representation  of  the  M2-F2  taken  from  Reference 
[48]  is  given  in  Figure  4-33.  The  upper  eleven  flaps  provided  coarse  longitudinal  trim  and  were 
deflected  differentially  (aileron  deflection)  to  provide  roll  control.  Because  of  adverse  yaw  due  to 
aileron  deflection,  a  mechanical  aileron-to-rudder  interconnect  was  provided.  The  interconnect 
ratio  was  adjustable  by  the  pilot.  A  stability  augmentation  system  was  used  to  improve  damping 
in  the  pitch,  roll,  and  yaw  axes. 

The  M2-F2  was  dropped  from  a  B-52  at  an  altitude  of 45,000  feet  and  a  Mach  number  of 
about  0.6.  Flight  test  maneuvers  were  performed  during  the  unpowered  descent  to  assess  the 
handling  characteristics  of  the  vehicle.  The  pilots  were  thoroughly  familiar  with  the  flight  plan 
and  predicted  handling  qualities  of  the  research  vehicle  as  a  result  of  training  on  a  six-degree-of- 
freedom  fixed  base  simulator.  The  M2-F2  exhibited  good  handling  qualities  longitudinally,  but 
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Figure  4-33.  Three-View  Drawing  of  the  M2-F2.  Dimensions  in  meters  (feet). 


the  lateral/directional  characteristics  were  rated  consistently  lower  than  the  longitudinal;  due 
partly  to  the  appearance  of  a  coupled  roll  spiral  mode  at  low  angles  of  attack.  The  vehicle 
experienced  severe  lateral/directional  oscillations  during  three  of  its  sixteen  flights.  Each  of  the 
oscillations  occurred  when  the  pilot  was  attempting  to  closely  control  the  bank  angle  at  negative 
angles  of  attack.  The  PIO  of  the  final  flight  contributed  to  a  gear-up  landing  in  which  the  vehicle 
was  extensively  damaged.  1^*1  As  the  result  of  the  accident,  a  center  fin  was  added  to  the  vehicle 
to  improve  the  lateral/directional  dynamics  and  the  slightly  redesigned  vehicle  was  designated  the 
M2-F3. 


4.3.7. 1  Develop  Pilot  Model 

The  task  under  consideration  is  the  tight  control  of  the  bank  angle  and  will  be  modeled  as  a 
compensatory  tracking  task.  A  modified  optimal  control  model  of  the  human  pilot  will  be 
considered  for  this  example.  The  pilot  model  is  designed  based  on  the  linear  transfer  function  of 
the  controlled  system.  For  simplification  and  comparison  to  Reference  [48],  only  the  pilots 
aileron  input  is  considered.  This  restriction  should  not  pose  a  problem  as  simulation  analysis  in 
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Reference  [48]  indicated  PIO  for  the  situation  of  aileron  only  input  as  well  as  for  coordinated 
aileron/rudder  input  during  turn  to  final  approach.  A  simplified  block  diagram  of  the  situation 
under  consideration  is  given  in  Figure  4-34.  The  MOCM  code  of  Reference  [23]  was  used  to 
compute  the  pilot  model.  The  parameters  used  in  the  MOCM  model  are  identical  to  those  given 
in  Table  4-1  with  the  disturbance  filter  given  by  13.3  l/(s2+2(0.7)(0.5)s+(0.52)).  The  values  are 
based  on  those  given  in  Reference  [49]  for  the  roll  axis. 
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Figure  4-34.  Block  Diagram  for  Single  Axis  Tracking  Task. 


4.3. 7.2  Identification  of  the  Dynamic  Elements  Suspected  of  Contributing  to  PIO 

A  first  step  in  the  closed-loop  analysis  is  to  identify  those  elements  of  the  aircraft  and  stability 
augmentation  system  that  are  suspected  to  be  contributing  factors  in  the  PIO  occurrence.  A  block 
diagram  of  the  M2-F2  lateral/directional  flight  control  system  is  given  in  Figure  4-35.  Since 
actuator  rate  limiting  is  known  to  contribute  to  PIO,  the  rate  limits  of  both  the  aileron  and  rudder 
actuators  will  be  considered  in  the  analysis.  Furthermore,  since  a  mechanical  aileron-to-rudder 
interconnect  is  needed  to  provide  acceptable  lateral/directional  characteristics,  it  also  poses  a  * 

potential  threat  in  terms  of  oscillation  tendencies.  This  threat  was  recognized  during  the  fixed 
based  simulations  as  depicted  in  Figure  4-36  which  is  taken  from  Reference  [48].  Command 
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Figure  4-35.  Block  Diagram  of  M2-F2  Lateral/Directional  Control  System. 


K|,  deg/deg 


Figure  4-36.  Simulator  Predicted  Regions  of  Lateral  Control  Problems  as  a  Function  of  Angle  of 
Attack  and  Aileron-to-Rudder  Interconnect  Ratio. 
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authority  limits  were  placed  on  the  pilot  input  as  well  as  both  the  roll  and  yaw  SAS  commands. 
These  command  limits,  together  with  the  actuator  rate  limits,  comprise  a  total  of  five  nonlinear 
elements  that  will  be  considered  as  possible  contributors  to  PIO.  The  rate  and  command  limits 
are  given  in  Table  4-9.  Along  with  the  aileron-to-rudder  interconnect  ratio,  variations  in  the  roll 
and  yaw  SAS  gains  will  be  considered  as  linear  elements  of  the  vehicle  suspected  of  contributing 
to  undesirable  oscillations. 


Table  4-9.  M2-F2  Actuator  Rate  and  Command  Authority  Limits. 


Upper  Flap 

Rudder 

Pilot  Aileron 

Roll 

Yaw 

Rate 

Rate 

Command 

SAS 

SAS 

(Aileron) 

Authority 

Authority 

Authority 

30  deg/sec 

22  deg/sec 

+10  deg 

±5  deg 

+4.2  deg 

In  addition  to  elements  of  the  aircraft  and  SAS,  various  aspects  of  the  pilot  are  also 
considered.  The  underlying  concept  of  the  closed-loop  analysis  is  that  a  PIO  resistant  aircraft  will 
remain  stable  in  the  face  of  variations  in  the  pilot  model.  At  this  point,  certain  aspects  of  the 
model  seem  reasonable  to  consider.  The  most  obvious  suspect  is  the  pilot  gain.  Increased  pilot 
gain  is  often  associated  with  tight  tracking  tasks  such  as  approach  and  landing. 

4.3. 7.3  Stability  Robustness  Analysis 

Those  elements  identified  as  possible  contributors  to  PIO  in  Section  4.3. 7. 2  will  now  be 
considered  in  a  stability  robustness  analysis.  The  linear  elements  will  all  be  considered  to  have  a 
multiplicative  uncertainty  as  shown  in  Figure  4-3.  The  linear  elements  are  the  pilot  gain  (Kpiiot), 
the  roll  SAS  gain  (Kp),  the  yaw  SAS  gain  (K,),  and  the  interconnect  ratio  (Ki).  The  nonlinear 
elements  (the  three  command  limits  and  the  two  actuator  rate  limits)  will  be  modeled  with  the 
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equivalent  gain  representation  of  Figure  4-4.  The  SIMULINK™  diagram  used  to  define  M(s)  for 
the  stability  robustness  analysis  is  given  in  Figure  4-37.  Associated  with  each  uncertainty  to  be 
considered  is  an  inport  outport  pair.  The  inport/outport  pairs  and  the  corresponding  dynamic 
elements  are  listed  in  Table  4-10.  When  the  litmod  command  is  issued,  a  state-space 
representation  of  M(s)  is  returned  having  q  inputs  and  q  outputs  (q  equals  the  total  number  of 
uncertainties  so  q  =  6  in  this  case).  If  it  is  desired  to  consider  only  a  subset  of  the  defined 
uncertainties,  the  one  need  only  to  extract  the  relevant  portions  of  the  state-space  model.  If,  for 
example,  only  the  combined  effect  of  an  uncertain  pilot  gain  in  conjuncture  with  rate  limiting  of 
the  aileron  actuator  is  to  be  considered,  then  only  the  first  and  third  inputs  and  outputs  of  the 
system  M(s)  are  needed. 


Table  4-10.  Inport/Outport  Pairs  for  Stability  Robustness  Analysis. 


Inport/Outport 

Pair 


Description  of 
Dynamic  Element 


1 

2 

3 

4 

5 

6 


Aileron  Rate  Limit 
Rudder  Rate  Limit 
Pilot  Gain 
Roll  SAS  Gain 
Yaw  SAS  Gain 
Interconnect  Ratio 


With  the  M-A  construction  complete,  the  structured  singular  value  can  be  used  to  calculate 
the  pertinent  stability  margins.  As  stated  previously,  when  bounds  are  known  to  exist  on  a  given 
variation,  the  corresponding  uncertainty  should  be  scaled  to  reflect  this  variation.  If,  for  example, 
the  roll  SAS  gain  is  known  to  vary  by  no  more  than  10%  from  the  nominal  value,  then  the  scaling 
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Figure  4-37.  SIMULINK^^  Diagram  for  Stability  Robustness  Analysis  of  M2-F2. 


gain  of  0. 1  should  be  applied  to  the  corresponding  inport/outport.  This  scaling  ensures  that  as  6 
varies  from  -1  to  1,  the  corresponding  SAS  gain  varies  from  Kp(l- .  1)  to  Kp(l  +.  1),  or  a  variation 
of +10%  from  the  nominal  gain.  However  since  no  such  information  is  given  for  the  M2-F2,  an 
alternate  approach  is  used.  With  all  the  scaling  gains  set  to  unity,  the  structured  singular  value  for 
real  parameter  variations  is  used  to  estimate  the  single  uncertainty  stability  margins  for  each  of  the 
uncertainties.  The  results  of  this  analysis  are  presented  in  Table  4-11.  A  1  or  -1  in  a  column 
represents  the  corresponding  “worst  case”  variation  is  either  positive  or  negative.  The  value 
indicated  by  min(l/|.ir)  represents  the  size  of  the  real  parameter  variation  necessary  before 
instability  is  possible.  The  frequency  at  which  the  minimum  of  the  structured  singular  value 
inverse  occurs  is  denoted  (0mi„.  Also  given  in  the  table  are  the  exact  uncertainties  necessary  to 
cause  instability  and  the  corresponding  frequencies.  The  estimated  stability  margins  are  between 
93%  and  98%  of  the  exact  values. 


Table  4-11.  Single  Uncertainty  Stability  Margins  for  the  M2-F2. 


Pilot 

Gain 

Aileron 

Rate 

Rudder 

Rate 

Roll 

SAS 

Yaw 

SAS 

Interconnec 

Ratio 

Minimum 

(1/Pr) 

Sexact 

®mm 

(rad/sec) 

©exact 

(rad/sec) 

1 

0 

0 

0 

0 

0 

.6088 

.6315 

5.55 

5.92 

0 

-1 

0 

0 

0 

0 

.8571 

.9032 

4.11 

2.93 

0 

0 

-1 

0 

0 

0 

.6848 

.7285 

2.65 

2.28 

0 

0 

0 

-1 

0 

0 

1.7883 

1.9050 

4.78 

4.31 

0 

0 

0 

0 

-1 

0 

1.0485 

1.0516 

2.75 

2.70 

0 

0 

0 

0 

0 

-1 

.1325 

.1357 

3.30 

3.44 

For  those  uncertainties  with  min(l/|.ir)  <  1,  the  scaling  gains  are  now  set  equal  to  the  single 
uncertainty  stability  margin  estimates.  For  example,  the  scaling  gain  on  the  pilot  gain  uncertainty 
inport  is  set  to  0.6088.  This  results  in  scaling  the  uncertainties  to  their  single  loop  equivalent 
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stability  margins  with  the  exception  of  the  uncertainties  in  the  roll  and  yaw  SAS  gains  which  are 
left  at  unity.  Now,  the  effect  of  varying  the  parameters  simultaneously  can  be  considered  without 
undue  scaling  bias.  The  results  of  the  scaled  simultaneous  variations  are  given  in  Table  4-12.  The 
stability  margins  given  in  this  table  now  represent  margins  that  are  subordinate  to  the  single  loop 
margins.  Consider  the  case  of  varying  pilot  gain  and  the  interconnect  ratio  simultaneously.  The 
margin  is  given  as  0.6586.  This  means  that  the  closed-loop  system  is  guaranteed  to  remain  stable 
as  long  as  both  parameters  remain  below  66%  of  the  their  corresponding  single  loop  instability 
values.  This  means,  assuming  all  other  uncertainties  are  zero,  as  long  as  the  pilot  gain  does  not 
increase  by  more  than  0.6586*0.6088=40%  and  the  interconnect  ratio  does  not  decrease  by  more 
than  0.6586*0.1325=9%,  the  closed-loop  system  will  remain  stable. 


Table  4-12.  Multiple  Uncertainty  Stability  Margins  for  the  M2-F2. 


Pilot 

Gain 

Aileron 

Rate 

Rudder 

Rate 

Roll 

SAS 

Yaw 

SAS 

Interconnect 

Ratio 

Minimum 

(1/Hr) 

®min 

(rad/sec) 

1 

-1 

0 

0 

0 

0 

.6894 

5.12 

1 

0 

-1 

0 

0 

0 

.8196 

6.30 

1 

0 

0 

-1 

0 

0 

.6298 

5.24 

1 

0 

0 

0 

-1 

0 

.6736 

5.81 

1 

0 

0 

0 

0 

-1 

.6586 

3.01 

1 

-1 

-1 

0 

0 

0 

.6566 

5.62 

1 

-1 

0 

-1 

0 

0 

.5235 

5.06 

1 

0 

-1 

0 

-1 

0 

.5506 

2.88 

1 

0 

0 

-1 

-1 

0 

.4893 

5.55 

1 

-1 

0 

0 

1 

0 

.5559 

4.51 

1 

0 

-1 

-1 

0 

0 

.5879 

5.62 

1 

-1 

0 

0 

0 

1 

.5539 

3.79 

1 

0 

-1 

0 

0 

1 

.5076 

3.30 

1 

0 

0 

-1 

0 

-1 

.4961 

5.36 

1 

0 

0 

0 

-1 

1 

.4525 

3.12 

1 

-1 

0 

-1 

0 

-1 

.4351 

5.24 

1 

-1 

-1 

0 

0 

1 

.4662 

3.34 

1 

0 

-1 

0 

-1 

1 

.3696 

3.05 

1 

-1 

-1 

-1 

-1 

1 

.3277 

3.08 
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It  is  obvious  that  stability  margins  calculated  for  multiple  uncertainties  are  greatly  affected  by 
the  scaling  gains.  Therefore,  when  bounds  on  the  model  parameter  variations  are  not  known,  this 
analysis  can  only  be  used  to  determine  relative  stability  due  to  various  combinations  of 
uncertainties.  With  this  reasoning,  the  results  of  Table  4-12  can  be  viewed  as  a  sort  of  sensitivity 
analysis. 

Based  upon  a  percentage  from  nominal,  from  the  single  uncertainty  margins  of  Table  4-11,  the 
interconnect  ratio  distinguishes  itself  as  the  parameter  to  which  stability  of  the  closed-loop  system 
is  by  far  the  most  sensitive.  The  pilot  gain  is  a  distant  second.  Since  we  are  primarily  interested 
in  examining  the  stability  of  the  system  with  respect  to  variation  in  the  pilot,  all  of  the  multiple 
uncertainty  margins  given  in  Table  4-12  include  a  variation  in  the  pilot  gain.  The  first  5  rows  of 
Table  4-12  show  the  results  for  varying  only  two  parameters.  Again,  stability  of  the  system  is 
strongly  affected  by  variations  in  the  pilot  gain  and  the  interconnect  ratio.  The  next  six  rows 
show  the  results  for  varying  three  parameters  simultaneously  while  keeping  the  interconnect  ratio 
at  its  nominal  value.  The  remaining  rows  show  the  results  for  varying  pilot  gain,  interconnect 
ratio  and  some  combinations  of  the  other  parameters. 

Of  interest  to  note  are  the  results  for  varying  pilot  gain,  aileron  rate  limit,  and  roll  SAS  gain 
(row  #7)  as  compared  with  varying  pilot  gain,  rudder  rate  limit,  and  yaw  SAS  gain  (row  #8)  for 
the  cases  of  constant  (row  #1  and  #8)  and  varying  (row  #16  and  #18)  interconnect  ratios.  When 
the  interconnect  ratio  is  held  constant,  system  stability  is  more  sensitive  to  variations  in  the  aileron 
rate  limit  and  roll  SAS  gain  (min(l/|ir)  =  0  52  as  opposed  to  0.55).  System  stability  is  more 
sensitive  to  variations  in  the  rudder  rate  limit  and  the  yaw  SAS  gain  when  the  interconnect  ratio  is 
allowed  to  vary  (min(l/|.ir)  =  0.37  compared  to  0.44).  Since  the  pilot  command  affects  the  rudder 
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only  through  the  interconnect  ratio,  this  switch  in  sensitivity  is  understandable.  Note  also,  the 
“worst  case”  direction  of  variation  of  the  interconnect  ratio  is  different  for  the  two  cases. 

Another  trend  that  can  be  found  in  Table  4-12  is  that  an  increase  in  the  interconnect  ratio 
causes  a  decrease  in  the  corresponding  frequency.  This  trend  is  easily  seen  in  comparing  the  first 
four  rows  to  rows  12  through  15.  A  last  trend  to  note  is  perhaps  the  most  obvious.  Namely,  the 
stability  margins  decrease  as  any  single  parameter  is  added  to  an  already  existing  combination. 
The  stability  margin  for  varying  pilot  gain,  aileron  rate  limit  ,and  roll  SAS  gain  (0.5235)  is  less 
than  the  stability  margin  for  varying  only  pilot  gain  and  aileron  rate  limit  (0.6894), or  only  pilot 
gain  and  roll  SAS  gain  (0.6298). 

4.3. 7.4  Multivariable  Describing  Function  Analysis 

The  stability  robustness  analysis  of  the  previous  section  produced  guaranteed  stability  bounds. 

That  is  to  say,  as  long  as  the  given  parameters  stay  within  the  bounds,  the  closed-loop  system  will 
remain  stable  and  will  not  limit  cycle.  However,  the  structured  singular  value  analysis  yields  little 
information  about  the  behavior  of  the  instability  when  the  bounds  are  breached.  Also,  since  the 
variations  of  the  equivalent  gains  representing  the  nonlinearities  are  input  amplitude  dependent,  it 
is  more  difficult  to  ascertain  if  these  bounds  will  be  breached  by  the  operating  pilot/aircraft 
system.  To  provide  further  insight  into  the  system  behavior,  a  describing  function  analysis  is 
performed  to  identify  the  existence  of  any  limit  cycles. 

A  limit  cycle  can  only  exist  when  either  the  pilot,  aircraft,  or  both  are  modeled  as  nonlinear 
systems.  For  the  M2-F2  example,  the  pilot  is  assumed  to  be  the  linear  pilot  given  by  the  MOCM 
while  the  aircraft  is  assumed  to  contain  the  5  isolated  saturation  nonlinearities  stated  in  Section 
4. 3. 7. 3.  Figure  4-38  shows  the  SIMULINK™  diagram  used  to  determine  M(s)  for  the  limit  cycle 
search.  When  the  linmod  command  is  used,  a  state  space  representation  of  M(s)  is  returned 
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Figure  4-38.  SIMULINK™  Diagram  for  Describing  Function  Analysis  of  M2-F2. 


having  5  inputs  and  5  outputs,  one  input/output  pair  for  each  of  the  nonlinear  elements  modeled. 


Table  4-13  lists  the  SIMULINK^“  inport/outport  pairs  and  the  corresponding  dynamic  elements. 
If  it  is  desired  to  analyze  only  a  subset  of  the  nonlinearities  as  modeled,  the  new  system  cannot  be 
obtained  by  simply  examining  certain  elements  of  the  original  system  as  was  done  for  the  stability 
robustness  analysis.  Separate  SIMULINK^^  diagrams  containing  only  the  particular  subset  of  the 
nonlinearities  should  be  used  to  extract  the  desired  state  space  representation  of  M(s).  An  entirely 
new  M(s)  needs  to  be  generated  because  the  nominal  system  is  not  obtained  by  setting  A=0  in  the 
describing  function  analysis.  Instead,  the  loop  is  opened  at  each  nonlinear  element  when  A=0. 


Table  4-13.  Inport/Outport  Pairs  for  Describing  Function  Analysis. 


Inport/Outport 

Pair 


Description  of 
Dynamic  Element 


1 

2 

3 

4 

5 


Aileron  Rate  Limit 
Rudder  Rate  Limit 
Pilot  Command  Limit 
Roll  SAS  Command  Limit 
Yaw  SAS  Command  Limit 


With  the  system  recast  in  the  standard  M-A  form,  the  numerical  methods  of  Section  4.3. 4. 3 
can  be  used  to  search  for  the  existence  of  limit  cycles.  The  limit  cycle  search  is  a  search  for  the 
inputs  to  the  nonlinear  elements  of  the  form  of  equation  (4-17)  that  solve  the  harmonic  balance 
equation  (4-19).  This  search  can  include  many  of  the  same  variations  as  the  stability  robustness 
analysis.  Presented  here  will  be  the  results  of  considering  all  five  nonlinearities  in  conjunction 
with  a  variation  in  pilot  gain  and  a  variation  in  the  interconnect  ratio.  The  pilot  gain  and 
interconnect  are  suspected  of  contributing  to  the  oscillation  tendencies  as  is  indicated  by  the 
results  of  the  stability  robustness  analysis  of  the  previous  section.  The  results  of  the  limit  cycle 
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search  are  given  in  a  set  of  three  tables.  Each  table  contains  limit  cycle  solutions  over  a  range  of 
pilot  gains  for  a  single  value  of  the  interconnect  ratio. 

In  addition  to  the  inputs  to  the  nonlinear  elements,  the  magnitudes  and  phases  of  other 
signals  in  the  loop  may  be  of  interest.  In  the  present  case,  the  bank  angle  response  (<t»)  is  of 
primary  concern  as  an  indication  of  the  severity  of  the  oscillation.  When  other  such  signals  are 
desired,  they  may  be  approximated  by  standard  block  diagram  manipulation.  Figure  4-39  shows 
the  relationship  between  the  bank  angle  and  the  inputs  to  the  actuator  rate  limiting  nonlinearities. 
With  the  input  signals. 


xjt)  =  a,sin(cot  +  ej 
x^(t)  =  a,sin((ot  +0,) 

The  steady  state  bank  angle  output  is  given  by, 

m  =  a  .|G,(j®)|sin(cot +  0,  +  ZG.(jco)) +  a,|G,(jfl>)|sin((ot +  0,  +ZG,(j®)) 

=  a^sin((ot  +  0  J 

where  G,  and  Gr  are  as  indicated  in  Figure  4-39  and  the  magnitudes  and  phases  are  to  be 
calculated  at  the  limit  cycle  solution.  Na  and  H  are  SIDF  representations  of  the  saturation 
nonlinearity.  The  bank  angle  output  phase  and  amplitude  can  then  be  given  as. 


0^  =  tan' 

fa. 

G. 

sin(0,  +ZG,)+a,|G,|sin(0,  +ZG,)^ 

u. 

G.|cos(0„  +ZG.)  +  a,  GJcos(0,  +ZG,)j 

(4-60) 

(4-61) 


(4-62) 


(4-63) 


a.|G.|sin(0.  +  ZG.)  +  a,|G,|  sin(0,  +ZG,) 
sin0^ 

^  a.|G.|cos(0.  +ZG.)  +  a,|G,|cos(0,  +ZG,) 
COS0^ 

The  bank  angle  results  as  calculated  above  are  also  included  in  the  tables. 
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Figure  4-39.  Steady  State  Bank  Angle  Output  Due  to  Sinusoidal  Inputs. 

Table  4-14  gives  the  results  for  the  nominal  interconnect  ratio  of  0.45.  The  table  shows  limit 
cycle  solutions  exist  for  the  given  range  of  pilot  gains.  These  solutions  are  characterized  by  a  high 
degree  of  aileron  saturation.  The  input  to  the  aileron  rate  limiter  has  an  amplitude  between  five 
and  nine  times  the  limit  value.  The  only  other  saturation  element  that  is  limited  is  the  pilot’s 
aileron  command.  The  commanded  input  is  between  one  to  two  times  the  command  limit.  The 
remaining  elements  are  not  saturated.  The  solutions  all  oscillate  at  a  frequency  of  approximately  4 
rad/sec.  The  limit  cycle  solutions  are  predicted  to  be  stable  for  pilot  gains  greater  than  or  equal  to 
1.5  times  the  nominal  gain  and  unstable  at  lower  gains.  For  the  nominal  value  of  the  interconnect 
ratio,  the  bank  angle  amplitude  of  the  predicted  limit  cycles  is  seen  to  be  on  the  order  of  two 
degrees. 

From  the  stability  robustness  analysis,  the  gain  margin  for  varying  only  pilot  gain  was  found  to 
be  1.6  It  was  found  that  the  stable  limit  cycles  predicted  in  Table  4-14  exhibit  soft  self  excitation 
for  pilot  gains  greater  than  this  margin  and  hard  self  excitation  for  pilot  gains  below  the  margin. 
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In  fact,  all  stable  limit  cycles  found  for  the  M2-F2  with  MOCM  pilot  models  seem  to  follow  this 


trend.  Namely,  soft  self  excitation  is  observed  when  the  parameters  are  set  such  that  the  stability 


margins  of  the  linear  system  representation  are  violated. 


The  effect  of  varying  interconnect  ratio  on  the  predicted  limit  cycle  can  be  seen  in  Tables  4-14 
through  4-16.  As  shown  in  Table  4-15,  lowering  the  interconnect  ratio  to  0.375  causes  the  bank 
angle  amplitude  and  the  frequency  of  the  oscillation  to  increase  slightly.  Also  noted  is  that  stable 
limit  cycle  solutions  are  predicted  for  lower  pilot  gains  than  for  the  nominal  interconnect  ratio.  As 
the  ratio  is  increased  to  0.525  (Table  4-16),  the  predicted  bank  angle  amplitude  increases  to 
approximately  10  degrees  and  the  frequency  decreases  to  under  3  rad/sec.  For  the  increased 
interconnect  ratio,  stable  limit  cycles  are  predicted  for  the  entire  range  of  pilot  gains  given. 


Table  4-14.  Limit  Cycle  Solution  for  Nominal  Interconnect  Ratio  (M=0.45). 


Kp 


w 

(rad/e) 


02 

(rad) 


03 

(rad) 


04 

(rad) 


05 

(rad) 


bank  bank 
Amp.  phase 


stable 

da/dAl 


2.00 

8.435 

3.985 

0.839 

3.141 

1.794  -1.172 

1.90 

8.295 

4.016 

0.840 

3.140 

1.681  -1.176 

1.80 

8,115 

4.056 

0.842 

3.138 

1.565  -1.182 

1.70 

7.870 

4.111 

0.844 

3.136 

1.444  -1.191 

1.60 

7.510 

4.194 

0.848 

3.132 

1.313  -1.205 

1.50 

6.845 

4.353 

0.853 

3.124 

1.153  -1.232 

o 

1—1 

5.866 

4.533 

0.856 

3.120 

1.033  -1.289 

1.30 

5.707 

4.442 

0.852 

3.135 

1.044  -1.318 

1,20 

5.517 

4.328 

0,850 

3.151 

1.058  -1.353 

1.10 

5.282 

4.183 

0.850 

3.168 

1.077  -1.393 

1.00 

4.980 

3.999 

0.853 

3.181 

1.103  -1.442 

0.450  -2.394 

0.134 

2.518 

2.850 

2.176 

-0.535 

0.448  -2.389 

0.134 

2.505 

2.815 

2.182 

-0.477 

0.445  -2.383 

0.134 

2.488 

2.772 

2.189 

-0.414 

0.442  -2.377 

0.134 

2.465 

2.714 

2.198 

-0.343 

0.437  -2.368 

0.133 

2.431 

2.631 

2.209 

-0.257 

0.428  -2.357 

0.131 

2.367 

2.481 

2.225 

-0.136 

0.427  -2.380 

0.130 

2.341 

2.385 

2.212 

0.006 

0.447  -2.435 

0.137 

2.437 

2.571 

2.165 

0.014 

0.469  -2.504 

0.148 

2.518 

2.790 

2.105 

0.022 

0.492  -2.593 

0.164 

2.571 

3.058 

2.027 

0.029 

0.520  -2.711 

0.181 

2.586 

3.395 

1.921 

0.033 

Table  4-15.  Limit  Cycle  Solutions  for  Decreased  Interconnect  Ratio  (KI=0.375). 


Kp 

A1 

ro 

A2 

92 

A3  03 

A4  04 

A5 

05 

bank 

bank 

stable 

( rad/s ) 

(rad) 

(rad) 

(rad) 

(rad) 

Amp. 

phase 

da/dAl 

2.00 

9.732 

4.560 

0.798 

3.099 

1.954  -1.083 

0.574  -2.166 

0.186 

1.931 

3.170 

2.422 

-0.363 

1.90 

9.615 

4.575 

0.798 

3.098 

1.844  -1.088 

0.572  -2.168 

0.185 

1.926 

3.150 

2.420 

-0.345 

CD 

o 

9.470 

4.593 

0.799 

3.096 

1.733  -1.095 

0.570  -2.171 

0.185 

1.920 

3.126 

2.418 

-0.324 

1.70 

9.288 

4.617 

0.799 

3.094 

1.620  -1.103 

0.567  -2.174 

0.184 

1.912 

3.095 

2.415 

-0.299 

1.60 

9.052 

4.648 

0.800 

3.091 

1.504  -1.113 

0.564  -2.179 

0.183 

1.902 

3.055 

2.411 

-0.268 

1.50 

8.730 

4.692 

0.800 

3.088 

1.384  -1.128 

0.559  -2.186 

0.181 

1.887 

2.998 

2.405 

-0.230 

1.40 

8.251 

4.761 

0.801 

3.082 

1.254  -1.150 

0.551  -2.197 

0.179 

1.865 

2.912 

2.396 

-0.178 

1.30 

7.329 

4.905 

0.802 

3.069 

1.096  -1.196 

0.534  -2.219 

0.174 

1.820 

2.738 

2.377 

-0.090 

1.20 

6.342 

5.022 

0.794 

3.069 

1.008  -1.264 

0.532  -2.273 

0.163 

1.842 

2.690 

2.334 

0.004 

1.10 

6.196 

4.961 

0.783 

3.091 

1.009  -1.294 

0.555  -2.320 

0.158 

1.956 

2.887 

2.301 

0.011 

1.00 

6.046 

4.889 

0.774 

3.116 

1.010  -1,326 

0.580  -2.373 

0.159 

2.088 

3.108 

2.263 

0.017 
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Table  4-16.  Limit  Cycle  Solutions  for  Increased  Interconnect  Ratio  (KI=0.525). 


Kp  A1  frt  A2  02  A3  03  A4  04  A5  05  bank  bank  stable 


( rad/fl ) 

(rad) 

(rad) 

(rad) 

(rad) 

Amp. 

phase  dcr/dAl 

2.00 

5.905 

2.691 

0.848 

3.100 

5.231  -1.612 

1.043  -3.433 

0.120 

4.350 

9.907 

1.070  -0.545 

1.90 

5.879 

2.694 

0.848 

3.100 

4.956  -1.613 

1.040  -3.432 

0.120 

4.347 

9.867 

1.071  -0.485 

1.80 

5.849 

2.697 

0.848 

3.101 

4,679  -1.614 

1.036  -3.432 

0.120 

4.343 

9.820 

1.072  -0.425 

1.70 

5.811 

2.701 

0.848 

3.101 

4.402  -1.615 

1.032  -3.431 

0.120 

4.338 

9.763 

1.073  -0.367 

1.60 

5.764 

2.706 

0.848 

3.101 

4.122  -1.616 

1.026  -3.430 

0.120 

4.332 

9.692 

1.074  -0.309 

1.50 

5.703 

2.712 

0.848 

3.102 

3.840  -1.618 

1.019  -3.429 

0.120 

4.325 

9.604 

1.076  -0.250 

1.40 

5.620 

2.720 

0.848 

3.102 

3.554  -1.621 

1.010  -3.428 

0.120 

4.315 

9.491 

1.078  -0.190 

1.30 

5.498 

2.731 

0.848 

3.103 

3.262  -1.624 

0.998  -3.426 

0.120 

4.302 

9.336 

1.080  -0.108 

1.20 

5.313 

2.746 

0.848 

3.104 

2.962  -1,629 

0.981  -3.423 

0.119 

4.283 

9.118 

1.004  -0.107 

1.10 

5.052 

2.769 

0.848 

3.105 

2.650  -1.636 

0.956  -3.418 

0.  119 

4.256 

0.816 

1.090  -0.105 

1.00 

4.660 

2.804 

0.847 

3.107 

2.321  -1.647 

0.920  -3.412 

0.118 

4.216 

8.373 

1.100  -0.099 

A  graphical  representation  of  the  limit  cycle  solutions  for  varying  interconnect  ratio  for  a  fixed 
pilot  gain  of  two  times  the  nominal  gain  is  given  in  Figures  4-40  and  4-41.  Figure  4-40  shows  the 
frequency  of  the  oscillation  decreases  slightly  as  the  interconnect  ratio  is  increased.  Recall,  this 
trend  was  also  the  case  in  the  stability  robustness  analysis.  The  bank  angle  amplitude  is  seen  to  be 
a  minimum  at  the  nominal  interconnect  ratio  of  0.45.  The  bank  angle  amplitude  increases  slightly 
as  the  interconnect  ratio  is  decreased  from  its  nominal  value.  As  the  interconnect  ratio  is 
increased  from  the  nominal  value,  the  bank  angle  amplitude  grows  rapidly. 

The  corresponding  input  amplitudes  to  the  nonlinear  elements  are  shown  in  Figure  4-41 .  The 
input  amplitudes  have  been  scaled  such  that  values  greater  than  unity  indicate  saturation  of  the 
limiters.  The  input  amplitude  to  the  yaw  SAS  command  limiter  and  the  input  amplitude  to  the 
rudder  rate  limiter  are  always  less  than  unity  indicating  these  elements  are  not  saturated  during  the 
limit  cycles.  The  input  amplitude  to  the  aileron  command  limiter  and  the  input  amplitude  to  the 
aileron  rate  limiter  are  always  greater  than  unity  indicating  these  elements  are  saturated  during  the 
limit  cycles.  The  input  amplitude  to  the  aileron  rate  limiter  generally  decreases  with  increasing 
interconnect  ratio.  The  input  amplitude  to  the  aileron  command  limiter  is  a  minimum  at  the 
nominal  interconnect  ratio  and  follows  the  same  trend  as  the  bank  angle  amplitude.  The  bank 
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Figure  4-40.  Predicted  Bank  Angle  Amplitude  and  Frequency. 


Figure  4-41.  Predicted  Input  Amplitudes  to  the  Saturation  Elements. 
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angle  amplitude  is  approximately  twice  the  scaled  value  of  the  amplitude  of  the  input  to  the 
aileron  command  limiter. 

All  of  the  predicted  limit  cycles  indicated  rate  limiting  of  the  aileron  actuator  as  well  as 
limiting  of  the  pilots  command  input.  The  only  other  limiting  observed  was  that  of  the  roll  SAS 
command  when  the  interconnect  ratio  was  increased  above  0.52.  This  analysis  shows  the 
importance  of  the  aileron  rate  limit  and  the  pilots  command  authority  limit  on  predicted  oscillation 
events.  From  this  analysis  it  would  seem  reasonable  the  nonlinearities  associated  with  the  rudder 
rate  limit  and  the  yaw  SAS  could  have  been  neglected  in  the  analysis. 

The  preceding  analysis  predicts  the  existence  of  limit  cycles  as  the  solutions  to  the  harmonic 
balance  equation.  Computer  simulation  of  the  closed-loop  pilot/vehicle  system  can  then  be  used 
to  verify  the  results.  Figure  4-42  shows  the  bank  angle  response  for  the  three  interconnect  ratios 
given  in  Tables  4-14  through  4-16  with  the  pilot  gain  in  each  case  set  to  twice  the  nominal  value 
for  that  case.  The  plot  generally  verifies  the  frequency  and  amplitude  data  given  in  the  tables. 
However,  at  the  larger  interconnect  ratio  the  simulation  results  in  a  bank  angle  amplitude  that  is 
twice  that  predicted.  Keeping  in  mind  that  the  describing  function  analysis  is  an  approximation, 
and  remembering  the  bank  angle  amplitude  begins  to  increase  sharply  at  large  interconnect  ratios, 
the  discrepancy  is  understandable.  Figure  4-43  shows  the  closed  curve  on  the  bank  angle-bank 
angle  rate  phase  plane  that  is  characteristic  of  a  limit  cycle. 

For  a  given  set  of  parameter  values,  there  may  be  more  than  one  limit  cycle  solution.  Table 
4-17  shows  another  set  of  solutions  for  varying  pilot  gain  while  maintaining  the  nominal 
interconnect  ratio.  These  solutions  are  characterized  by  large  bank  angle  amplitudes  (80  degrees) 
and  low  frequencies  (1.2  rad/sec).  These  values  agree  much  better  with  the  flight  data  which 
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indicated  bank  angle  amplitudes  greater  than  50  degrees  at  a  frequency  of  approximately  1.6 
rad/sec.  The  six-degree-of-freedom  simulation  data  of  Reference  [48]  for  aileron  only  pilot 
control  indicated  bank  angle  amplitudes  near  80  degrees  at  a  frequency  of  1.3  rad/sec.  While 
these  values  agree  well  with  those  in  Table  4-17,  a  computer  simulation  of  this  predicted  limit 
cycle  has  yet  to  be  obtained.  Again  it  appears  that  the  phenomena  of  hard  self-excitation  is  the 
culprit. 


Table  4-17.  Additional  Limit  Cycle  Solutions  for  Nominal  Interconnect  Ratio  (KI=0.45). 


Kp 

A1 

A2 

02 

A3  03 

A4 

04 

A5  05 

bank 

bank 

stable 

( rad/ B ) 

(rad) 

(rad) 

(rad) 

(rad) 

Amp. 

phase 

da/dAl 

2.00 

0.627 

1.203 

0.391 

2.884 

12.513  -1.938 

3.584 

2.500 

0.384  -0.045 

82.427 

0.486 

-2.373 

l.BQ 

0.627 

1.203 

0.391 

2.884 

11.260  -1.938 

3.583 

2.500 

0.384  -0.045 

82.419 

0.486 

-2.379 

1.60 

0.627 

1.203 

0.391 

2.884 

10.007  -1.938 

3.583 

2.500 

0.384  -0.045 

82.407 

0.485 

-2.386 

1.40 

0.626 

1.203 

0.390 

2.884 

8.753  -1.938 

3.582 

2.500 

0.384  -0.045 

82.389 

0.485 

-2.396 

1.20 

0.626 

1.203 

0.390 

2.884 

7.499  -1.938 

3.580 

2.499 

0.384  -0.045 

82.362 

0.485 

-2.413 

1.00 

0.625 

1.203 

0.390 

2.884 

6.244  -1.939 

3.577 

2.499 

0.384  -0.045 

82.317 

0.485 

-2.441 

O.BO 

0.624 

1.202 

0.389 

2.883 

4.988  -1.940 

3.573 

2.499 

0.384  -0.045 

82.234 

0.484 

-2.495 

0.60 

0.621 

1.202 

0.387 

2.883 

3.729  -1.942 

3.562 

2.497 

0.383  -0.046 

82.051 

0.483 

-2.626 

0.40 

0.612 

1.199 

0.381 

2.881 

2.461  -1.947 

3.530 

2.494 

0.381  -0.047 

81.504 

0.478 

-3.144 

0.20 

0.542 

1.181 

0.338 

2.867 

1.138  -1.997 

3.290 

2.463 

0.364  -0.060 

77.394 

0.442 

2.266 

0.00 

0.181 

0.887 

0.088 

2.581 

0.000  -1.997 

2.086 

1.593 

0.339  -0.668 

69.958 

-0.550 

12.056 

4.3.7. 5  Summary 

The  stability  robustness  analysis  of  the  M2-F2  confirmed  that  the  interconnect  ratio,  aileron 
rate  limit,  and  pilot  gain  had  a  considerable  affect  on  the  closed-loop  system  stability.  (The 
analysis  also  indicates  rudder  rate  limiting  has  an  impact,  but  the  analysis  considers  no  pilot  rudder 


input,  so  it  is  unlikely  that  the  rudder  rate  limit  will  be  reached.)  In  particular,  closed-loop 
stability  is  very  sensitive  to  the  value  of  the  interconnect  ratio. 


The  describing  function  analysis  shows  that  self  sustaining  oscillations  of  the  pilot/aircraft 


system  are  possible  for  parameter  variations  in  the  range  of  interest.  Furthermore,  the 


dependency  of  the  character  of  these  oscillations  on  the  parameter  variations  emerges  through  the 


limit  cycle  analysis.  Though  the  predicted  characteristics  of  the  oscillations  that  can  currently  be 
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simulated  do  not  agree  well  with  that  experienced  in  flight,  this  result  is  not  cause  for  alarm. 

First,  the  linear  characteristics  of  the  M2-F2  are  so  poor  that  nonlinear  analysis  is  not  necessary  to 
predict  handling  qualities  problems,  including  PIO.  Further,  the  information  available  for  the  M2- 
F2  did  not  include  actuator  models  for  aileron  and  rudder.  Also,  the  analysis  was  carried  out  at  a 
constant  angle  of  attack  since  no  longitudinal  dynamics  were  given. 

One  very  important  point  needs  to  be  stated.  The  current  analysis  is  heavily  dependent  on  the 
pilot  model  used.  More  precisely,  the  analysis  is  dependent  on  the  pilot  model  and  parameter 
variations  of  that  model.  At  this  point  it  seems  logical  that  variations  in  pilot  gain  and  time  delay 
will  figure  prominently  in  future  analysis.  Though  some  idea  exists  about  what  to  vary,  more 
research  needs  to  be  performed  to  determine  how  much  of  a  variation  to  consider.  Since 
acceptable  ranges  of  the  parameter  variations  have  not  yet  been  established,  it  would  be  premature 
at  this  point  to  make  any  concrete  assertion  as  to  PIO  proneness  of  any  single  vehicle  based  on 
the  current  analysis. 


K 
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5  Conclusions  and  Recommendations 


The  new  unified  PIO  analysis  theory  shows  significant  promise  and  warrants  continued 
development.  It  can  handle  many  simultaneous  dynamic  effects  in  a  systematic  way  It  has  ties  to 
existing  criteria  and  it  can  be  used  to  isolate  combinations  of  individual  dynamic  effects  that  are 
causing  the  PIO. 

The  analysis  of  HAVE  PIO  configurations  revealed  that  the  PIO  tendencies  of  these 
configurations  can  be  predicted  using  simple  pilot  gain/phase  variation  arguments.  The  vector 
margin  concept,  drawn  from  multivariable  stability  theory,  correctly  identified  all  but  one  of  the 
HAVE  PIO  configurations.  These  configurations  were  also  used  to  draw  connections  between 
several  existing  PIO  criteria  and  the  new  unified  method. 

The  M2-F2  example  demonstrated  how  many  individual  dynamic  effects  can  be  analyzed 
simultaneously.  In  this  case,  the  combined  effect  of  five  separate  nonlinearities  and  two  linear 
parameter  variations  were  studied.  While  the  PIO  tendency  in  this  case  was  most  likely  caused  by 
poor  flying  qualities,  a  unified  PIO  analysis  shows  exactly  which  dynamic  characteristics  couple 
with  the  pilot  to  cause  oscillations.  For  example,  the  analysis  results  showed  that  the  PIO 
amplitude  and  frequency  that  occurred  during  flight  testing  and  manned  simulation  of  the  M2-F2 
could  be  accurately  predicted  by  assuming  an  increased  pilot  gain  and  simultaneous  limiting  in  the 
pilot  command  and  roll  SAS  feedback  channels. 

5.1  Plan  for  Follow  on  Work 

There  are  several  activities  that  should  be  carried  out  before  the  new  unified  PIO  analysis 
method  will  be  completely  accepted  in  industry.  Some  of  these  activities  are  intended  to  make  the 
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analysis  methods  easier  to  use  and  more  reliable.  There  are  a  few  remaining  research  questions  to 
address  as  well. 

One  remaining  research  question  relates  to  the  study  of  which  pilot  model  parameters  should 
be  varied  and  to  what  degree  they  should  be  varied.  Several  variations  in  pilot  gain,  time  delay, 
and  lead/lag  time  constants  were  studied  during  this  research.  While  pilot  gain  variations  seemed 
to  provide  the  best  PIO  predictions,  there  is  sufficient  evidence  to  warrant  further  investigation  of 
other  parameters.  For  example,  the  single-loop  analysis  of  the  HAVE  PIO  configurations  showed 
that  the  critical  frequency  did  not  exactly  predict  the  flight  test  PIO  frequency.  It  is  likely  that  a 
better  correlation  in  frequency  could  be  obtained  if  a  specific  variation  in  pilot  time  delay  were 
used. 

The  other  remaining  research  question  focuses  on  the  existence  of  hard  and  soft  self¬ 
excitations.  Recall  that  a  soft  self-excitation  results  when  a  sustained  limit  cycle  can  be  obtained 
by  simulating  an  initial  condition  response  only.  A  hard  self-excitation  requires  application  and 
removal  of  an  external  input  to  yield  a  sustained  oscillation.  During  this  research,  none  of  the 
predicted  hard  self-excitations  were  confirmed  using  the  SIMULINK  simulation  program.  Both 
random  and  deterministic  input  signals  were  tried.  Continued  investigation  of  the  hard  self¬ 
excitation  phenomenon  is  important  because  it  might  shed  light  on  the  fact  that  some  PIOs  occur 
even  after  the  aircraft  has  been  in  the  fleet  for  several  years. 

Successful  application  to  additional  case  studies  will  also  aid  in  industry  acceptance  of  the  new 

•t 

unified  method.  This  report  demonstrated  application  to  the  HAVE  PIO  configurations,  the  F-4 
*  (with  back-up  control  systems),  the  M2-F2,  and  one  flight  condition  of  the  B-2.  Analysis  of  other 

known  PlO-free  and  PlO-prone  aircraft  is  needed  to  establish  more  definite  design  guidelines. 
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Finally,  the  software  developed  for  this  work  is  of  research  quality  only.  There  are  many 
improvements  that  could  be  made  in  order  to  make  the  software  easier  to  use  and  perhaps  more 
general  in  application.  There  are  also  several  recently  published  algorithms  that  will  improve  the 
accuracy  of  the  structured  singular  value  calculations  used  in  the  unified  PIO  theory,  t^^l  Reliable 
and  efficient  software  developed  specifically  for  unified  PIO  analysis  will  undoubtedly  improve 
industry  acceptance. 
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